options(warn = -1)
options(repr.plot.width=15, repr.plot.height=15)
library(readxl)
col_I_sn<-read_excel("stat2021_sm\\us_col.std\\I_shortname.xls")
## New names:
## • `` -> `...1`
col_I<-read_excel("stat2021_sm\\us_col.std\\I.xls")
## New names:
## • `` -> `...1`
head(col_I_sn)|>rmarkdown::paged_table()
summary(col_I_sn)
## ...1 PPIND FICE STATE
## Length:176 Min. :1.000 Min. : 1009 Length:176
## Class :character 1st Qu.:1.000 1st Qu.: 1738 Class :character
## Mode :character Median :1.000 Median : 2550 Mode :character
## Mean :1.301 Mean : 2901
## 3rd Qu.:2.000 3rd Qu.: 3417
## Max. :2.000 Max. :10366
##
## TYPE AVRMATH AVRVERB AVRCOMB
## Length:176 Min. :390.0 Min. :391.0 Min. : 810
## Class :character 1st Qu.:520.0 1st Qu.:454.5 1st Qu.: 980
## Mode :character Median :544.0 Median :479.5 Median :1017
## Mean :563.5 Mean :494.7 Mean :1058
## 3rd Qu.:603.0 3rd Qu.:518.5 3rd Qu.:1126
## Max. :750.0 Max. :665.0 Max. :1410
## NA's :68 NA's :68 NA's :67
## AVR_ACT MATH_1 MATH_3 VERB_1
## Min. :19.00 Min. :350.0 Min. :480.0 Min. :320.0
## 1st Qu.:22.00 1st Qu.:460.0 1st Qu.:590.0 1st Qu.:400.0
## Median :23.00 Median :500.0 Median :627.5 Median :430.0
## Mean :23.69 Mean :515.4 Mean :635.0 Mean :448.3
## 3rd Qu.:25.00 3rd Qu.:570.0 3rd Qu.:687.5 3rd Qu.:480.0
## Max. :31.00 Max. :740.0 Max. :780.0 Max. :630.0
## NA's :83 NA's :34 NA's :34 NA's :34
## VERB_3 ACT_1 ACT_3 APP_REC
## Min. :440.0 Min. :16.00 Min. :21.00 Min. : 787
## 1st Qu.:520.0 1st Qu.:19.00 1st Qu.:25.00 1st Qu.: 4323
## Median :550.0 Median :21.00 Median :26.00 Median : 7654
## Mean :561.3 Mean :21.28 Mean :26.55 Mean : 8544
## 3rd Qu.:600.0 3rd Qu.:23.00 3rd Qu.:28.00 3rd Qu.:11776
## Max. :720.0 Max. :29.00 Max. :33.00 Max. :48094
## NA's :34 NA's :67 NA's :67 NA's :1
## APP_ACC NEW_STUD NEW10 NEW25 FULLTIME
## Min. : 507 Min. : 210 Min. : 8.00 Min. :24.00 Min. : 912
## 1st Qu.: 3033 1st Qu.:1264 1st Qu.:24.00 1st Qu.:52.00 1st Qu.: 5846
## Median : 4761 Median :1949 Median :32.00 Median :63.00 Median :10215
## Mean : 5546 Mean :2252 Mean :41.48 Mean :65.55 Mean :11296
## 3rd Qu.: 7232 3rd Qu.:3035 3rd Qu.:57.00 3rd Qu.:82.00 3rd Qu.:15033
## Max. :26330 Max. :7425 Max. :98.00 Max. :99.00 Max. :31643
## NA's :1 NA's :16 NA's :26
## PARTTIME IN_STATE OUT_STAT R_B_COST
## Min. : 16.0 Min. : 647 Min. : 2279 Min. :2082
## 1st Qu.: 804.8 1st Qu.: 2100 1st Qu.: 6712 1st Qu.:3588
## Median : 1694.5 Median : 3030 Median : 8400 Median :4213
## Mean : 2519.9 Mean : 6641 Mean :10108 Mean :4538
## 3rd Qu.: 3240.8 3rd Qu.:12348 3rd Qu.:12668 3rd Qu.:5564
## Max. :21836.0 Max. :20100 Max. :20100 Max. :7425
## NA's :6 NA's :7 NA's :1
## ROOM BOARD ADD_FEE BOOK PERSONAL
## Min. :1033 Min. :1000 Min. : 20.0 Min. : 300.0 Min. : 300
## 1st Qu.:1810 1st Qu.:1762 1st Qu.: 210.0 1st Qu.: 500.0 1st Qu.:1164
## Median :2644 Median :2125 Median : 425.5 Median : 600.0 Median :1600
## Mean :2708 Mean :2241 Mean : 648.1 Mean : 603.1 Mean :1763
## 3rd Qu.:3490 3rd Qu.:2586 3rd Qu.: 694.0 3rd Qu.: 673.8 3rd Qu.:2200
## Max. :6965 Max. :4760 Max. :4374.0 Max. :1230.0 Max. :6800
## NA's :42 NA's :64 NA's :40 NA's :2 NA's :11
## PH_D TERM_D SF_RATIO DONATE
## Min. :63.00 Min. :67.00 Min. : 2.90 Min. : 4.00
## 1st Qu.:80.50 1st Qu.:87.00 1st Qu.:10.88 1st Qu.:10.75
## Median :87.00 Median :92.00 Median :14.50 Median :17.00
## Mean :85.72 Mean :90.52 Mean :14.23 Mean :19.01
## 3rd Qu.:92.00 3rd Qu.:96.00 3rd Qu.:18.02 3rd Qu.:24.00
## Max. :99.00 Max. :99.00 Max. :24.70 Max. :54.00
## NA's :9 NA's :16 NA's :12
## INSTRUCT GRADUAT SAL_FULL SAL_AC
## Min. : 3605 Min. :10.00 Min. : 446.0 Min. :364.0
## 1st Qu.: 7604 1st Qu.:47.50 1st Qu.: 597.0 1st Qu.:445.0
## Median : 9840 Median :62.00 Median : 661.0 Median :479.5
## Mean :12832 Mean :62.02 Mean : 669.2 Mean :487.1
## 3rd Qu.:14340 3rd Qu.:74.50 3rd Qu.: 732.2 3rd Qu.:521.2
## Max. :62469 Max. :99.00 Max. :1009.0 Max. :733.0
## NA's :5
## SAL_AS SAL_ALL COMP_FUL COMP_AC
## Min. :323.0 Min. :362.0 Min. : 537.0 Min. :438.0
## 1st Qu.:381.8 1st Qu.:472.2 1st Qu.: 729.0 1st Qu.:556.5
## Median :407.0 Median :522.5 Median : 815.0 Median :606.0
## Mean :412.8 Mean :534.0 Mean : 827.2 Mean :612.1
## 3rd Qu.:434.2 3rd Qu.:578.2 3rd Qu.: 910.0 3rd Qu.:662.2
## Max. :576.0 Max. :866.0 Max. :1236.0 Max. :909.0
##
## COMP_AS COMP_ALL NUM_FULL NUM_AC
## Min. :395.0 Min. : 436.0 Min. : 39.0 Min. : 32.0
## 1st Qu.:481.0 1st Qu.: 587.0 1st Qu.:184.5 1st Qu.:138.5
## Median :508.5 Median : 652.0 Median :278.0 Median :208.0
## Mean :519.0 Mean : 665.8 Mean :336.8 Mean :230.3
## 3rd Qu.:552.2 3rd Qu.: 730.2 3rd Qu.:457.8 3rd Qu.:299.0
## Max. :717.0 Max. :1075.0 Max. :997.0 Max. :721.0
## NA's :1
## NUM_AS NUM_INS NUM_ALL
## Min. : 29.0 Min. : 0.00 Min. : 108.0
## 1st Qu.:128.5 1st Qu.: 5.00 1st Qu.: 505.5
## Median :175.0 Median : 16.00 Median : 721.0
## Mean :190.7 Mean : 27.09 Mean : 812.8
## 3rd Qu.:238.2 3rd Qu.: 40.00 3rd Qu.:1035.0
## Max. :510.0 Max. :178.00 Max. :2261.0
## NA's :1
col_I_sn$PPIND=ifelse(col_I_sn$PPIND==1, "public", "private")
NA review
col_I_sn_noNA<-col_I_sn|>na.omit()
paste("col_I_sn dim: ", paste(dim(col_I_sn),collapse=", "))
## [1] "col_I_sn dim: 176, 49"
paste("col_I_sn_noNA dim: ", paste(dim(col_I_sn_noNA), collapse=", "))
## [1] "col_I_sn_noNA dim: 22, 49"
naCol<-colSums(is.na(col_I_sn))|>sort(decreasing=TRUE)
naCol
## AVR_ACT AVRMATH AVRVERB AVRCOMB ACT_1 ACT_3 BOARD ROOM
## 83 68 68 67 67 67 64 42
## ADD_FEE MATH_1 MATH_3 VERB_1 VERB_3 NEW25 NEW10 TERM_D
## 40 34 34 34 34 26 16 16
## DONATE PERSONAL PH_D IN_STATE PARTTIME GRADUAT BOOK APP_REC
## 12 11 9 7 6 5 2 1
## APP_ACC OUT_STAT NUM_AC NUM_INS ...1 PPIND FICE STATE
## 1 1 1 1 0 0 0 0
## TYPE NEW_STUD FULLTIME R_B_COST SF_RATIO INSTRUCT SAL_FULL SAL_AC
## 0 0 0 0 0 0 0 0
## SAL_AS SAL_ALL COMP_FUL COMP_AC COMP_AS COMP_ALL NUM_FULL NUM_AS
## 0 0 0 0 0 0 0 0
## NUM_ALL
## 0
table(naCol)
## naCol
## 0 1 2 5 6 7 9 11 12 16 26 34 40 42 64 67 68 83
## 21 5 1 1 1 1 1 1 1 2 1 4 1 1 1 3 2 1
log dataset
q_cols<-sapply(col_I_sn, is.numeric)
q_cols[c("PPIND", "FICE")]<-FALSE
q_cols<-names(q_cols)[q_cols==TRUE]
col_I_sn_log<-col_I_sn|>
mutate(across(q_cols, ~log(.x)))|>
rename_with(~paste0("log-", .x), q_cols)
head(col_I_sn_log)|>rmarkdown::paged_table()
Выделим основные признаки, по которым будет проходить дальнейший анализ данных.
Всего 49 признаков.
Отберём из условий задачи и по смыслу
Из условий задачи: 2. ADD_FEE - дополнительные выплаты 3. BOOK - плата за бронирование 4. NEW10 (зависимая переменная, также NEW25)
NEW10 зависит от AVRMATH, MATH_1, MATH_3, AVRVERB, VERB_1, VERB_3, AVR_ACT, ACT_1, ACT_3 и AVRCOMB. Признаков много, так что выберем самый обобщающий - AVRCOMB 5. AVRCOMB
SF_RATIO
PH_D
GRADUAT 6-8 из предположения, что они взаимосвязаны
R_B_COST
IN_STATE
OUT_STAT 9-11 - зависимость стоимости проживания (как окажется дальше из matrix plot, IN_STATE - фиктивный признак, то есть не информативный)
Проделаем это для рассматриваемых 11 признаков.
Setup
colNames<-c("ADD_FEE", "BOOK", "R_B_COST","IN_STATE","OUT_STAT","NEW10","AVRCOMB","SF_RATIO","PH_D","GRADUAT","PPIND")
colFin<-c("ADD_FEE", "BOOK", "PPIND")
colCost<-c("R_B_COST", "IN_STATE", "OUT_STAT", "PPIND")
colNew<-c("NEW10", "AVRCOMB", "PPIND")
colStud<-c("SF_RATIO", "PH_D", "GRADUAT", "PPIND")
col_I_comp<-col_I_sn[colNames]
col_I_comp_log<-col_I_sn_log[c(paste0("log-", colNames)[-length(colNames)], "PPIND")]
head(col_I_comp_log)|>rmarkdown::paged_table()
unique_info <- function(column) {
column<-na.omit(column)
list(ratio=length(unique(column)) / length(column), unique=length(unique(column)), total=length(column), moda = max(table(column)))
}
ratios <- sapply(col_I_comp, unique_info)
ratios
## ADD_FEE BOOK R_B_COST IN_STATE OUT_STAT NEW10 AVRCOMB
## ratio 0.875 0.3793103 0.9431818 0.9053254 0.9314286 0.425 0.7889908
## unique 119 66 166 153 163 68 86
## total 136 174 176 169 175 160 109
## moda 4 29 3 4 4 10 3
## SF_RATIO PH_D GRADUAT PPIND
## ratio 0.6363636 0.1976048 0.3918129 0.01136364
## unique 112 33 67 2
## total 176 167 171 176
## moda 5 10 7 123
ggpairs(col_I_comp, aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
col_I_sn|>filter(ADD_FEE>3000)
University of California, Massachusetts - особые индивиды, в связи с высоким ВВП штатов.
Видим, что NEW10 зависит от PH_D, GRADUAT, AVRCOMB, отдельно для private зависит от OUT_STAT, IN_STATE, для public - от ADD_FEE, BOOK, R_B_COST
Прологарифмируем ADD_FEE
col_I_mixed<-col_I_comp|>mutate(ADD_FEE=col_I_comp_log$`log-ADD_FEE`)|>rename(`log-ADD_FEE`=ADD_FEE)|>dplyr::select(-IN_STATE)
head(col_I_mixed)
ggpairs(col_I_mixed, aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Особые индивиды
col_I_sn|>filter(BOOK>=900 | BOOK <= 350)
Howard, Tulane, Hofstra University, Rensselaer Polytechn - особенности расположения, штата.
ggpairs(col_I_mixed[,c("log-ADD_FEE", "BOOK", "PPIND")], aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
col_I_sn|>filter(BOOK>=900 | BOOK <= 350)
Может быть как низкой, так и высокой в зависимости от штата, политики университета и скрытых факторов.
col_I_sn|>filter(ADD_FEE>3000)
University of California, Massachusetts - особые индивиды (из-за ВВП).
ggpairs(col_I_mixed[colNew], aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
col_I_sn|>filter(AVRCOMB<=1200 & NEW10>=75)
University of California, North Carolina. Причиной может быть низкий порог для поступления.
ggpairs(col_I_mixed[colStud], aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
col_I_sn|>filter(PH_D<=70)
Незначительные выбросы - особые индивиды. Могут быть связаны с низким GRADUAT и высоким SF_RATIO.
ggpairs(col_I_mixed[,c("R_B_COST", "OUT_STAT", "PPIND")], aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
col_I_sn|>filter(OUT_STAT>=14000 & PPIND=="public")
col_I_sn|>filter(OUT_STAT<=10500 & PPIND=="private")
Цены на обучение варьируются по штатам и самим университетам.
ggpairs(col_I_mixed, aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Выделим неоднородности R_B_COST, OUT_STAT. Явное разделение на два облака точек для этих признаков может быть связано с разницей финансирования. Для public - государство, штат и граждане, для private - только граждане.
results <- describeBy(col_I_mixed, group=col_I_comp$PPIND)
df_results <- map_dfr(results, ~as.data.frame(.x), .id = "group")
kable(df_results, "html") %>%
kable_styling(full_width = FALSE) %>%
scroll_box(width = "100%", height = "400px")
| group | vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| log-ADD_FEE…1 | private | 1 | 46 | 5.838785 | 0.7936613 | 6.003811 | 5.872662 | 0.7083206 | 3.688880 | 7.515345 | 3.826465 | -0.4493357 | 0.0057261 | 0.1170190 |
| BOOK…2 | private | 2 | 53 | 610.075472 | 158.1109456 | 600.000000 | 597.534884 | 148.2600000 | 300.000000 | 1230.000000 | 930.000000 | 1.2325851 | 3.0965764 | 21.7182087 |
| R_B_COST…3 | private | 3 | 53 | 5839.811321 | 997.0704793 | 5975.000000 | 5922.255814 | 816.9126000 | 3320.000000 | 7425.000000 | 4105.000000 | -0.7233444 | -0.2295138 | 136.9581633 |
| OUT_STAT…4 | private | 4 | 53 | 15747.358491 | 3924.9053349 | 17020.000000 | 16306.186047 | 2772.4620000 | 2340.000000 | 20100.000000 | 17760.000000 | -1.3082365 | 1.3192370 | 539.1272103 |
| NEW10…5 | private | 5 | 52 | 57.788461 | 24.3796255 | 58.000000 | 57.761905 | 33.3585000 | 16.000000 | 98.000000 | 82.000000 | 0.0381849 | -1.4295240 | 3.3808458 |
| AVRCOMB…6 | private | 6 | 29 | 1177.034483 | 144.2860459 | 1199.000000 | 1177.680000 | 161.6034000 | 883.000000 | 1410.000000 | 527.000000 | -0.0893581 | -1.2296542 | 26.7932461 |
| SF_RATIO…7 | private | 7 | 53 | 9.813208 | 4.3268854 | 9.200000 | 9.555814 | 4.8925800 | 2.900000 | 20.500000 | 17.600000 | 0.4758614 | -0.7282486 | 0.5943434 |
| PH_D…8 | private | 8 | 49 | 89.591837 | 7.8658502 | 91.000000 | 90.317073 | 7.4130000 | 71.000000 | 99.000000 | 28.000000 | -0.8078233 | -0.4898133 | 1.1236929 |
| GRADUAT…9 | private | 9 | 52 | 78.711539 | 15.7534142 | 78.500000 | 80.023809 | 17.7912000 | 33.000000 | 99.000000 | 66.000000 | -0.6127581 | -0.1790876 | 2.1846055 |
| PPIND…10 | private | 10 | 53 | 1.000000 | 0.0000000 | 1.000000 | 1.000000 | 0.0000000 | 1.000000 | 1.000000 | 0.000000 | NaN | NaN | 0.0000000 |
| log-ADD_FEE…11 | public | 1 | 90 | 6.014273 | 1.1296149 | 6.093567 | 6.024473 | 0.8723448 | 2.995732 | 8.383433 | 5.387701 | -0.1980167 | 0.2896088 | 0.1190719 |
| BOOK…12 | public | 2 | 121 | 600.008264 | 101.1621220 | 600.000000 | 591.103093 | 111.1950000 | 400.000000 | 858.000000 | 458.000000 | 0.6447669 | -0.2889320 | 9.1965565 |
| R_B_COST…13 | public | 3 | 123 | 3977.073171 | 842.4904744 | 3811.000000 | 3911.979798 | 722.0262000 | 2082.000000 | 6607.000000 | 4525.000000 | 0.7886999 | 0.7902446 | 75.9648078 |
| OUT_STAT…14 | public | 4 | 122 | 7657.549180 | 2252.3124723 | 7446.500000 | 7517.163265 | 1939.2408000 | 2279.000000 | 15732.000000 | 13453.000000 | 0.7771996 | 1.2996173 | 203.9147900 |
| NEW10…15 | public | 5 | 108 | 33.620370 | 21.8669373 | 26.000000 | 29.920455 | 11.8608000 | 8.000000 | 95.000000 | 87.000000 | 1.5420644 | 1.5799521 | 2.1041470 |
| AVRCOMB…16 | public | 6 | 80 | 1014.662500 | 84.6586510 | 997.000000 | 1010.687500 | 75.6126000 | 810.000000 | 1240.000000 | 430.000000 | 0.4180047 | 0.1997793 | 9.4651249 |
| SF_RATIO…17 | public | 7 | 123 | 16.129268 | 3.8302021 | 16.500000 | 16.178788 | 4.1512800 | 6.700000 | 24.700000 | 18.000000 | -0.0994206 | -0.5699091 | 0.3453577 |
| PH_D…18 | public | 8 | 118 | 84.110169 | 7.7985139 | 85.000000 | 84.406250 | 7.4130000 | 63.000000 | 99.000000 | 36.000000 | -0.3698537 | -0.3490528 | 0.7179114 |
| GRADUAT…19 | public | 9 | 119 | 54.731092 | 15.5730632 | 54.000000 | 54.412371 | 16.3086000 | 10.000000 | 95.000000 | 85.000000 | 0.1127780 | -0.1892703 | 1.4275804 |
| PPIND…20 | public | 10 | 123 | 2.000000 | 0.0000000 | 2.000000 | 2.000000 | 0.0000000 | 2.000000 | 2.000000 | 0.000000 | NaN | NaN | 0.0000000 |
NEW10 и AVRCOMB имеют skew, близкий к 0 (private), SF_RATIO, GRADUAT, log-ADD_FEE околонулевой skew (public).
Близкий к 0 (private): log-ADD_FEE (но skew ~ -0.44).
col_split_mixed<-split(col_I_mixed, col_I_mixed$PPIND)
names(col_split_mixed)<-c("private", "public")
colMixed<-names(col_I_mixed)
colMixed<-colMixed[-length(colMixed)]
colMixed
## [1] "log-ADD_FEE" "BOOK" "R_B_COST" "OUT_STAT" "NEW10"
## [6] "AVRCOMB" "SF_RATIO" "PH_D" "GRADUAT"
chi_squared_normality_test <- function(data) {
data <- na.omit(data)
n <- length(data)
if (n < 10) stop("Sample size too small for this test. Need at least 10 observations.")
mu <- mean(data)
sigma <- sd(data)
# Calculate quantiles for breaks
probs <- seq(5/n, 1 - 5/n, by = 5/n)
breaks <- qnorm(probs, mean = mu, sd = sigma)
breaks <- c(qnorm(0), breaks, qnorm(1))
#Observed Frequencies
observed <- hist(data, breaks = breaks, plot = FALSE)$counts
#Expected Frequencies (using pnorm)
expected <- diff(pnorm(breaks, mean = mu, sd = sigma)) * n
chi_squared <- sum((observed - expected)^2 / expected)
df <- length(expected) - 3 #Corrected degrees of freedom
p_value <- pchisq(chi_squared, df = df, lower.tail = FALSE)
return(list(p_value = p_value, statistic = chi_squared, df = df))
}
n = 500
test<-replicate(5000, chi_squared_normality_test(rnorm(n))$p_value)
quantile(test)
## 0% 25% 50% 75% 100%
## 0.000116153 0.249925966 0.509619748 0.761571168 0.999679814
perform_normality_tests <- function(df, cols) {
# Input validation (remains the same)
# Function to perform a single normality test (simplified)
run_test <- function(data, test_name) {
if (test_name == "Chi-squared") {
chisq_test<-chi_squared_normality_test(data)
return(chisq_test[c("p_value", "statistic", "df")])
}
test_result <- switch(test_name,
"Lilliefors" = lillie.test(data),
"Anderson-Darling" = ad.test(data),
"Shapiro-Wilk" = shapiro.test(data)
)
return(list(p_value = test_result$p.value, statistic = test_result$statistic, df=length(na.omit(data)))) # Removed test_name
}
results <- lapply(cols, function(col) {
data <- df[[col]]
n_na <- sum(is.na(data))
if (n_na > 0) {
warning(paste0("Removed ", n_na, " NA values from column '", col, "' before testing."))
data <- na.omit(data)
}
test_results <- lapply(c("Lilliefors", "Anderson-Darling", "Shapiro-Wilk", "Chi-squared"), function(test) run_test(data, test))
list(column = col, tests = test_results)
})
results_df <- do.call(rbind, lapply(results, function(x) {
test_names <- c("Lilliefors", "Anderson-Darling", "Shapiro-Wilk", "Chi-squared")
p_values <- sapply(x$tests, "[[", "p_value")
data.frame(
Column = x$column,
Test = test_names, #Directly use test_names vector
Statistic = sapply(x$tests, "[[", "statistic"),
P_value = p_values,
Size = sapply(x$tests, "[[", "df"),
Significance = sapply(p_values, function(p) ifelse(p < 0.05, "\u2714", "\u2716")),
stringsAsFactors = FALSE
)
}))
plots <- lapply(cols, function(col) {
data <- df[[col]]
data <- na.omit(data)
ggqqplot(data, title = paste("Normal Probability Plot for", col))
})
names(plots)<-cols
list(results = results_df, plots = plots)
}
# colSums(is.na(col_split_mixed$public))
col_norm_public<-col_split_mixed$public|>perform_normality_tests(colMixed)
col_norm_public$plots
## $`log-ADD_FEE`
##
## $BOOK
##
## $R_B_COST
##
## $OUT_STAT
##
## $NEW10
##
## $AVRCOMB
##
## $SF_RATIO
##
## $PH_D
##
## $GRADUAT
col_norm_public$results|>dplyr::select(-Statistic)|>rmarkdown::paged_table()
Хи-квадрат: не отвергаем AVRCOMB, SF_RATIO, PH_D, GRADUAT, R_B_COST, OUT_STAT,
Можно заметить “хвосты” на некоторых графиках.
col_norm_private<-col_split_mixed$private|>perform_normality_tests(colMixed)
col_norm_private$plots
## $`log-ADD_FEE`
##
## $BOOK
##
## $R_B_COST
##
## $OUT_STAT
##
## $NEW10
##
## $AVRCOMB
##
## $SF_RATIO
##
## $PH_D
##
## $GRADUAT
col_norm_private$results|>dplyr::select(-Statistic)|>rmarkdown::paged_table()
Хи-квадрат: не отвергаем log-ADD_FEE, R_B_COST, NEW10, AVRCOMB, SF_RATIO, GRADUAT
Совместно со следующим заданием.
perform_tests <- function(data,
quantitative_columns,
categorical_column,
alpha = 0.05) {
results <- list()
for (col in quantitative_columns) {
levene_test<-leveneTest(data[[col]] ~ data[[categorical_column]])
var_test <- var.test(data[[col]] ~ data[[categorical_column]])
wilcox_test <- wilcox.test(data[[col]]~data[[categorical_column]])
t_test_equal <- t.test(data[[col]] ~ data[[categorical_column]], var.equal = TRUE)
t_test_unequal <- t.test(data[[col]] ~ data[[categorical_column]], var.equal = FALSE)
ks_test<-ks.test(data[[col]]~data[[categorical_column]])
results[[col]] <- list(
t_test_unequal_p_value = t_test_unequal$p.value,
t_test_equal_p_value = t_test_equal$p.value,
wilcox_test_p_value = wilcox_test$p.value,
levene_test_p_value = levene_test$`Pr(>F)`[1],
var_test_p_value = var_test$p.value,
ks_test_p_value = ks_test$p.value
)
}
return(results |> enframe(name = "attribute", value = "p_values") |>
unnest_wider(p_values))
}
df_tests = perform_tests(col_I_mixed, colMixed, "PPIND")
df_long <- col_I_mixed|>
pivot_longer(cols = -PPIND, names_to = "attribute", values_to = "value")
# Объединяем результаты теста с основными данными
df_with_tests <- df_long %>%
left_join(df_tests, by = "attribute")
df_with_tests_cost <- df_with_tests[df_with_tests$attribute%in%colCost[-length(colCost)],]
df_with_tests_new <- df_with_tests[df_with_tests$attribute%in%colNew[-length(colNew)],]
df_with_tests_stud <- df_with_tests[df_with_tests$attribute%in%colStud[-length(colStud)],]
df_with_tests_fin <- df_with_tests[df_with_tests$attribute%in%c("log-ADD_FEE", "BOOK"),]
p<-ggplot(df_with_tests_cost,
aes(x = attribute,
y = value,
fill = PPIND)) +
geom_boxplot() +
labs(title = "Boxplot",
x = "Attribute",
y = "Value")+
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), plot.caption = element_text(size = 8, hjust = 0)) +
scale_fill_manual(values = c("private" = "red", "public" = "lightblue"))
p$data<-df_with_tests_new|>filter(attribute%in%"AVRCOMB")
p
p$data<-df_with_tests_new|>filter(attribute%in%"NEW10")
p
df_with_tests|>dplyr::filter(attribute%in%colNew[-length(colNew)])|>dplyr::select(-PPIND, -value)|>slice_head(n=2)|>rmarkdown::paged_table()
Можем сказать, что public и private для AVRCOMB и NEW10 имеют разные матожидания.
В случае гипотезы о равенстве дисперсий, критерий Левена (Levene) и критерий Фишера (var) имеют разные мощности. Критерий Фишера используется, если сравниваемые выборки имеют нормальное распределение. Поскольку группа public признака NEW10 отвергла гипотезу о нормальном распределении, применять критерий Фишера нельзя.
Отвергаем равенство дисперсий для AVRCOMB, NEW10.
p$data<-df_with_tests_fin|>filter(attribute%in%"BOOK")
p
p$data<-df_with_tests_fin|>filter(attribute%in%"log-ADD_FEE")
p
df_with_tests|>dplyr::filter(attribute%in%c("log-ADD_FEE", "BOOK"))|>dplyr::select(-PPIND, -value)|>head(length(colFin)-1)|>rmarkdown::paged_table()
Исходя из p_value не отвергаем, что public и private имеют равные матожидания. Тест Вилкоксона также характеризует их однородность.
Проверим log-ADD_FEE и BOOK друг с другом
compare_column_tests <- function(column_1,
column_2,
alpha = 0.05) {
levene_test<-leveneTest(column_2, column_1)
var_test <- var.test(column_1, column_2)
wilcox_test <- wilcox.test(column_1, column_2)
t_test_equal <- t.test(column_1, column_2, var.equal = TRUE)
t_test_unequal <- t.test(column_1, column_2, var.equal = FALSE)
results <- data.frame(
t_test_unequal_p_value = t_test_unequal$p.value,
t_test_equal_p_value = t_test_equal$p.value,
levene_test_p_value = levene_test$`Pr(>F)`[1],
var_test_p_value = var_test$p.value,
wilcox_test_p_value = wilcox_test$p.value
)
}
af_pr<-col_I_sn$ADD_FEE[col_I_sn$PPIND=="private"]
af_pu<-col_I_sn$ADD_FEE[col_I_sn$PPIND=="public"]
b_pr<-col_split_mixed$private$BOOK
b_pu<-col_split_mixed$public$BOOK
print(compare_column_tests(af_pr, b_pr))
## t_test_unequal_p_value t_test_equal_p_value levene_test_p_value
## 1 0.005903872 0.003538701 2.925758e-75
## var_test_p_value wilcox_test_p_value
## 1 1.42501e-07 1.678377e-05
print(compare_column_tests(af_pu, b_pu))
## t_test_unequal_p_value t_test_equal_p_value levene_test_p_value
## 1 0.1706573 0.1115642 0.9495775
## var_test_p_value wilcox_test_p_value
## 1 0 0.0006091651
В общем случае тесты показывают, что разница есть в случае общественных вузов, но для частных вузов не отвергаем равенство дисперсий (из levene’s test) и матожиданий.
p$data<-df_with_tests_stud
p
df_with_tests|>dplyr::filter(attribute%in%colStud[-length(colStud)])|>dplyr::select(-PPIND, -value)|>head(length(colStud)-1)|>rmarkdown::paged_table()
Отвергаем равенство матожиданий для GRADUAT, SF_RATIO, PH_D. Не отвергаем равенство дисперсий GRADUAT, SF_RATIO, PH_D. — ### OUT_STAT, R_B_COST
p$data<-df_with_tests_cost
p
df_with_tests|>dplyr::filter(attribute%in%colCost[-length(colCost)])|>dplyr::select(-PPIND, -value)|>head(length(colCost)-2)|>rmarkdown::paged_table()
Отвергаем равенство матожиданий для OUT_STAT, R_B_COST. Не отвергаем равенство дисперсий для R_B_COST.
t-критерий неудовлетворителен в нескольких ситуациях:
В качестве замены t-критерия можно воспользоваться непараметрическим критерием Манна-Уитни (Вилкоксона).
df_tests|>dplyr::select(attribute,t_test_unequal_p_value, t_test_equal_p_value, wilcox_test_p_value)|>rmarkdown::paged_table()
В случаях, где Манн-Уитни мощнее t-test: BOOK, R_B_COST, OUT_STAT, SF_RATIO, GRADUAT; причиной этого могут быть выбросы. В случаях, где Манн-Уитни не мощнее t-test: log-ADD_FEE, NEW10, AVRCOMB, PH_D; могут отсутствовать выбросы, из-за чего t-test оказывается мощнее.
Как непараметрический критерий, тест Вилкоксона уступает t-тесту в мощности против альтернативы о неравенстве матожиданий. С другой стороны, тест Вилкоксона устойчив к выбросам, мощный против альтернативы о неравенстве распределения двух случайных величин. Если известно, что две случайные величины имеют одинаковые формы распределения, то тест Вилкоксона автоматически ставит гипотезу о равенстве матожиданий.
df_tests|>dplyr::select(attribute, wilcox_test_p_value, ks_test_p_value)|>rmarkdown::paged_table()
Столбец wilcox_test_p_value - p-value из критерия Вилкоксона. Столбец ks_test_p_value - p-value из критерия Колмогорова-Смирнова. Критерии Вилкоксона и Колмогорова-Смирнова мощны против альтернативы о неравенстве распределения. Критерии не отвергают равенство форм распределений для log-ADD_FEE и BOOK, но отвергают для остальных. Среди рассматриваемых признаков нет случаев, когда один критерий отвергает гипотезу, а другой - не отвергает.
col_split_mixed<-split(col_I_mixed, col_I_mixed$PPIND)
names(col_split_mixed)<-c("private", "public")
col_split_mixed<-lapply(col_split_mixed, function(x){
dplyr::select(x, -PPIND)
})
colMixed<-names(col_I_mixed)
colMixed<-colMixed[-length(colMixed)]
col_split_mixed$private|>rmarkdown::paged_table()
col_split_mixed$public|>rmarkdown::paged_table()
ggpairs(col_split_mixed$private, aes(alpha = 0.5),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
Можно заметить, что NEW10/AVRCOMB разбивается на две части (поскольку AVRCOMB имеет много пропусков и сам AVRCOMB подразумевает разделение на подгруппы).
col_I_sn[col_I_sn$PPIND=="private",]|>filter(NEW10<=50 & !is.na(AVRCOMB))|>rmarkdown::paged_table()
col_I_sn[col_I_sn$PPIND=="private",]|>filter(NEW10>50 & !is.na(AVRCOMB))|>rmarkdown::paged_table()
Среди частных университетов только один “Delaware University” имеет расхождения между IN_STATE и OUT_STAT. Сейчас идёт такое ценообразование.
University of Delaware’s tuition is $15,410 for in-state and $37,930 for out-of-state.
col_I_sn[col_I_sn$IN_STATE!=col_I_sn$OUT_STAT & col_I_sn$PPIND=="private",]
ggpairs(col_split_mixed$public, aes(alpha = 0.5),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
На pairs plot можно заметить выбросы в признаках log-ADD_FEE, OUT_STAT, NEW10, которые влияют на коэффициент корреляции Пирсона.
Учитывая ограниченный объём данных, тем более после разделения по группам (например, наименьшая в private размера 29 без NA), следует воспользоваться pairwise. Более того, из-за признака AVRCOMB использование casewise неявно рассматривает корреляцию в другой подгруппе. Минусом pairwise является то, что из-за попарной проверки могут учитываться данные, пропущенные в других парах.
compute_corr <- function(group_data) {
corr_pairwise_pearson <- cor(group_data, use = "pairwise.complete.obs", method = "pearson")
corr_pairwise_spearman <- cor(group_data, use = "pairwise.complete.obs", method = "spearman")
return(list(pairwise_pearson = corr_pairwise_pearson, pairwise_spearman = corr_pairwise_spearman))
}
make_corr_plots <- function(corr_list, group) {
p1 <- ggcorrplot(corr_list$pairwise_pearson,
lab = TRUE,
lab_size = 1.9,
colors = c("blue", "white", "red"),
title = paste("Pairwise Cor Pearson Matrix for Group", group),
ggtheme = theme_minimal() +
theme(plot.title = element_text(size = 6)) # Reduced title size
)
p2 <- ggcorrplot(corr_list$pairwise_spearman,
lab = TRUE,
lab_size = 1.9,
colors = c("blue", "white", "red"),
title = paste("Pairwise Cor Spearman Matrix for Group", group),
ggtheme = theme_minimal() +
theme(plot.title = element_text(size = 6)) # Reduced title size
)
grid.arrange(p1, p2, ncol = 2, widths = c(4, 4), # Adjust widths here
padding = unit(1, "cm")) # Add padding for extra space
}
library(ggcorrplot)
for (group in names(col_split_mixed)) {
group_data <- col_split_mixed[[group]]
cor_matrix <- compute_corr(group_data)
plot_pearson_spearman <- make_corr_plots(cor_matrix, group)
print(plot_pearson_spearman)
}
## TableGrob (1 x 2) "arrange": 2 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## TableGrob (1 x 2) "arrange": 2 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
Пирсон мощнее Спирмена:
Спирмен мощнее Пирсона:
Пирсон мощнее Спирмена:
Спирмен мощнее Пирсона:
В тех ячейках, где совпадают или близки коэффициенты корреляции Пирсона и Спирмена почти нет выбросов и у scatterplot эллиптическая структура, если Пирсон значительно больше Спирмена - есть заметные выбросы. Спирмен может быть больше Пирсона, если выбросы имеются, но их положение не очень значимо относительно других наблюдений.
Значительные корреляции:
Задача: провести регрессию качества поступающих студентов (NEW10) на остальные признаки. С помощью пошаговой регрессии уменьшить число признаков
library(lm.beta)
Рассматриваемые наблюдения: 1. PPIND (для группировки) 2. ADD_FEE - дополнительные выплаты 3. BOOK - плата за бронирование 4. NEW10 (зависимая переменная, также NEW25) - студентов из топ 10% своей старшей школы 5. AVRCOMB - средний комбинированный балл 6. SF_RATIO - Student/Faculty ratio 7. PH_D - Pct. of faculty with Ph.D.’s 8. GRADUAT - Graduation rate 9. R_B_COST - Room and board costs 10. OUT_STAT - Out-of-state tuition
colNames<-c("...1", "ADD_FEE", "BOOK", "R_B_COST","OUT_STAT","NEW10","AVRCOMB","SF_RATIO","PH_D","GRADUAT","PPIND")
col_I_regr<-col_I_sn[c(colNames)]
col_I_regr<-col_I_regr|>mutate(ADD_FEE = log(ADD_FEE))|>rename(log_ADD_FEE=ADD_FEE)
col_I_regr_split<-split(col_I_regr, col_I_regr$PPIND)
# Автор: -Николай-
library(GGally)
library(ggplot2)
library(dplyr)
library(rlang)
##
## Присоединяю пакет: 'rlang'
## Следующие объекты скрыты от 'package:purrr':
##
## %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, splice
plot_numeric_pairs <- function(df, title = NULL, highlight = NULL) {
# Выбираем имена всех numeric-столбцов
numeric_cols <- df %>% dplyr::select(dplyr::where(is.numeric)) %>% names()
# Функция для нижних панелей с возможностью выделения точек разными цветами
custom_points <- function(data, mapping, ...) {
# Отрисовываем базовые точки синим цветом
p <- ggplot(data, mapping) +
geom_point(alpha = 0.6, size = 2, color = "blue", ...)
if (!is.null(highlight)) {
# Если highlight не является списком, превращаем его в список с одним элементом и цветом "red"
if (!is.list(highlight) || inherits(highlight, "data.frame")) {
highlight <- list(red = highlight)
}
# Получаем имена переменных из mapping с помощью as_label()
xvar <- as_label(mapping$x)
yvar <- as_label(mapping$y)
# Для каждого набора точек в списке highlight
for (col in names(highlight)) {
hl_df <- highlight[[col]]
# Проверяем, что hl_df содержит нужные переменные
if(all(c(xvar, yvar) %in% names(hl_df))) {
# Находим совпадающие строки между данными панели и hl_df
data_highlight <- semi_join(data, hl_df, by = c(xvar, yvar))
# Накладываем выделенные точки нужного цвета
p <- p + geom_point(data = data_highlight, mapping = mapping,
alpha = 0.6, size = 2, color = col, ...)
}
}
}
p
}
# Строим матрицу парных графиков
ggpairs(df,
columns = numeric_cols,
columnLabels = numeric_cols,
title = title,
upper = list(continuous = wrap("cor", size = 3)),
lower = list(continuous = custom_points),
diag = list(continuous = wrap("barDiag", bins = 10, fill = "lightblue", color = "black"))
) +
theme_bw() +
theme(
strip.text.x = element_text(size = 7, angle = 0, hjust = 0),
strip.text.y = element_text(size = 8, angle = 0, hjust = 0),
axis.text = element_text(size = 6),
plot.margin = unit(c(1, 2, 1, 1), "cm")
)
}
pr_data<-col_I_regr_split$private
print(paste0("Число наблюдений: ", dim(pr_data)[1]))
## [1] "Число наблюдений: 53"
colSums(is.na(pr_data))
## ...1 log_ADD_FEE BOOK R_B_COST OUT_STAT NEW10
## 0 7 0 0 0 1
## AVRCOMB SF_RATIO PH_D GRADUAT PPIND
## 24 0 4 1 0
Много пропусков в AVRCOMB
plot_numeric_pairs(pr_data, "Частные вузы")
outlier_pr_book_eye<-pr_data[pr_data$BOOK>875,]
outlier_pr_out_eye<-pr_data[pr_data$OUT_STAT<5000,]
outlier_pr_book_eye
outlier_pr_out_eye
plot_numeric_pairs(pr_data, "Частные вузы", highlight = list(red=outlier_pr_book_eye, black=outlier_pr_out_eye))
Здесь условия из .tsk файла
Стоит помнить, что у AVRCOMB много NA. В такой модели не будет выбросов “на глаз”
model_pr_new <- lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT,
data = pr_data)|>lm.beta()
summary(model_pr_new)
##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST +
## PH_D + SF_RATIO + GRADUAT + OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.4148 -3.8409 0.0271 2.9131 12.4176
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -78.090327 NA 51.790742 -1.508 0.1575
## AVRCOMB 0.117625 0.544789 0.042417 2.773 0.0169 *
## BOOK -0.006699 -0.029617 0.014003 -0.478 0.6410
## log_ADD_FEE 3.200735 0.098598 2.219759 1.442 0.1749
## R_B_COST -0.003437 -0.117762 0.002499 -1.376 0.1941
## PH_D -0.492675 -0.126378 0.496416 -0.992 0.3406
## SF_RATIO -1.353492 -0.205756 0.594526 -2.277 0.0419 *
## GRADUAT 0.306142 0.178421 0.258574 1.184 0.2594
## OUT_STAT 0.002354 0.257544 0.001391 1.693 0.1163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.797 on 12 degrees of freedom
## (32 пропущенных наблюдений удалены)
## Multiple R-squared: 0.9646, Adjusted R-squared: 0.9411
## F-statistic: 40.92 on 8 and 12 DF, p-value: 1.496e-07
# Автор: -Вадим-
create_ellipse <- function(center, cov_matrix, level = 0.95, n = 100) {
# Уровень доверия преобразуем в квантиль хи-квадрат
chi_sq <- qchisq(level, df = 2)
# Вычисляем собственные значения и собственные векторы
eigen_res <- eigen(cov_matrix)
# Углы для точек эллипса
angles <- seq(0, 2 * pi, length.out = n)
# Радиусы эллипса (основаны на собственных значениях)
radii <- sqrt(chi_sq * eigen_res$values)
# Создаем эллипс в стандартной ориентации
ellipse_points <- cbind(radii[1] * cos(angles), radii[2] * sin(angles))
# Поворачиваем эллипс в соответствии с собственными векторами
ellipse_points <- t(eigen_res$vectors %*% t(ellipse_points))
# Добавляем смещение
ellipse_points <- sweep(ellipse_points, 2, center, "+")
ellipse_points <- as.data.frame(ellipse_points)
colnames(ellipse_points) <- c("x", "y") # Assign column names
return(ellipse_points)
}
plot_confidence_ellipse <- function(model, which.coeff = c(2, 4), level = 0.95,
xlab = NULL, ylab = NULL, main = "Доверительный эллипс",
margin = 0.1) {
# Проверка на то, что переданы именно 2 коэффициента
if (length(which.coeff) != 2) {
stop("which.coeff должен содержать ровно 2 индекса коэффициентов")
}
# Получение необходимых данных из модели
model_summary <- summary(model)
coef_estimates <- coef(model)
сov_matrix <- model_summary$cov.unscaled
# Проверка допустимости индексов коэффициентов
if (any(which.coeff > length(coef_estimates))) {
stop("Индекс коэффициента превышает их количество в модели")
}
# Создание эллипса
ellipse_points <- create_ellipse(
center = coef_estimates[which.coeff],
cov_matrix = сov_matrix[which.coeff, which.coeff],
level = level
)
# Определение границ графика с учетом отступов
x_min <- min(ellipse_points$x)
x_max <- max(ellipse_points$x)
y_min <- min(ellipse_points$y)
y_max <- max(ellipse_points$y)
x_range <- c(x_min - margin * (x_max - x_min), x_max + margin * (x_max - x_min))
y_range <- c(y_min - margin * (y_max - y_min), y_max + margin * (y_max - y_min))
# Если не указаны названия осей, используем имена коэффициентов
if (is.null(xlab)) {
xlab <- names(coef_estimates)[which.coeff[1]]
}
if (is.null(ylab)) {
ylab <- names(coef_estimates)[which.coeff[2]]
}
# Создаем data.frame с центральной точкой
center_point <- data.frame(
x = coef_estimates[which.coeff[1]],
y = coef_estimates[which.coeff[2]]
)
# Создаем объект ggplot
p <- ggplot2::ggplot() +
# Добавляем эллипс
ggplot2::geom_path(data = ellipse_points, ggplot2::aes(x = x, y = y), color = "blue") +
# Добавляем точку с оценкой коэффициента
ggplot2::geom_point(data = center_point, ggplot2::aes(x = x, y = y), color = "red", size = 3) +
# Добавляем вертикальную и горизонтальную линии через центр
ggplot2::geom_vline(xintercept = center_point$x, linetype = "dashed", color = "gray") +
ggplot2::geom_hline(yintercept = center_point$y, linetype = "dashed", color = "gray") +
# Настраиваем оси и заголовок
ggplot2::xlim(x_range) +
ggplot2::ylim(y_range) +
ggplot2::labs(
x = xlab,
y = ylab,
title = main
) +
ggplot2::theme_light()
return(p)
}
plot_confidence_ellipse(model_pr_new, c(3, 2), level=0.95)
На примере этой пары: AVRCOMB и SF_RATIO значимы
На примере, PH_D и OUT_STAT
phcor<-cor(pr_data$PH_D, pr_data$OUT_STAT, use = "complete.obs")
y <- pr_data$NEW10
X <- pr_data[, c("PH_D", "OUT_STAT")]
model_cor <- lm(NEW10 ~ PH_D + OUT_STAT, data = pr_data)
sum_cor<-summary(model_cor)
sum_cor
##
## Call:
## lm(formula = NEW10 ~ PH_D + OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.696 -13.720 0.979 15.266 34.815
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.009e+02 3.392e+01 -2.975 0.00470 **
## PH_D 1.551e+00 4.763e-01 3.257 0.00214 **
## OUT_STAT 1.198e-03 9.551e-04 1.255 0.21608
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.96 on 45 degrees of freedom
## (5 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4428, Adjusted R-squared: 0.4181
## F-statistic: 17.88 on 2 and 45 DF, p-value: 1.926e-06
# Оцениваем дисперсию ошибок
sigma_s_squared <- (sum_cor$sigma/sd(y, na.rm=TRUE))^2 # оценка дисперсии ошибок
# Размер выборки
n <- nobs(model_cor)
R_XX<-cor(X, use="complete.obs")
print("Inverse R_XX")
## [1] "Inverse R_XX"
solve(R_XX)
## PH_D OUT_STAT
## PH_D 2.088331 -1.507579
## OUT_STAT -1.507579 2.088331
cov_beta<-sigma_s_squared/n*solve(R_XX)
print("Cov_beta:")
## [1] "Cov_beta:"
cov_beta
## PH_D OUT_STAT
## PH_D 0.02360939 -0.01704377
## OUT_STAT -0.01704377 0.02360939
print("Determinant")
## [1] "Determinant"
detcoef<-1/(1-phcor^2)
detcoef
## [1] 2.088331
print("-rho/(1-rho^2)")
## [1] "-rho/(1-rho^2)"
-phcor*detcoef
## [1] -1.507579
phcor
## [1] 0.7219063
# Автор: -Евгений-
# Вычисление множественной корреляции
calc_multicollinearity <- function(var, data) {
X <- as.matrix(data[, setdiff(names(data), var)]) # Матрица остальных регрессоров
y <- as.matrix(data[[var]]) # Текущий регрессор как зависимая переменная
model <- lm(y ~ X)
sqrt(summary(model)$r.squared)
}
calc_partial_correlation <- function(var, data) {
# Остатки регрессии текущего регрессора на остальные регрессоры
data<-data|>na.omit()
model_X <- lm(data[[var]] ~ ., data = data[, !names(data) %in% c(var, "NEW10")])
res_X <- residuals(model_X)
# Остатки регрессии Y на остальные регрессоры
model_Y <- lm(NEW10 ~ ., data = data[, !names(data) %in% var])
res_Y <- residuals(model_Y)
cor(res_X, res_Y, use = "complete.obs") # Корреляция остатков
}
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)
## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.6803712 0.3842868
## BOOK BOOK 0.4955889 -0.1368050
## R_B_COST R_B_COST 0.8079355 -0.3690491
## OUT_STAT OUT_STAT 0.9472357 0.4390291
## AVRCOMB AVRCOMB 0.9764603 0.6249356
## SF_RATIO SF_RATIO 0.8649231 -0.5492085
## PH_D PH_D 0.9121799 -0.2754193
## GRADUAT GRADUAT 0.9401219 0.3234121
Сильно положительные множественные корреляции вместе с частными корреляциями, следует сравнить по отдельности
regr_high <- c("R_B_COST", "OUT_STAT", "AVRCOMB", "SF_RATIO", "PH_D", "GRADUAT")
regr_fullhigh<-c(regr_high, "NEW10")
redund_high <- data.frame(
Регрессор = regr_high,
Множественная_корреляция = sapply(regr_high, calc_multicollinearity, pr_data[, regr_fullhigh]),
Частная_корреляция = sapply(regr_high, calc_partial_correlation, pr_data[,regr_fullhigh])
)
print(redund_high)
## Регрессор Множественная_корреляция Частная_корреляция
## R_B_COST R_B_COST 0.6621454 -0.1634685
## OUT_STAT OUT_STAT 0.9214250 0.3142483
## AVRCOMB AVRCOMB 0.9687653 0.5706967
## SF_RATIO SF_RATIO 0.8369497 -0.3982266
## PH_D PH_D 0.8830917 -0.1158525
## GRADUAT GRADUAT 0.9332980 0.3403546
Видно, что PH_D имеет высокую множественную корреляцию и низкую частную корреляцию, поэтому выбросим её
model_pr_noPhD<-lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr_noPhD)
##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST +
## SF_RATIO + GRADUAT + OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.9908 -2.5824 0.6409 3.1699 11.8038
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.156e+02 NA 2.833e+01 -4.080 0.000872 ***
## AVRCOMB 1.291e-01 6.067e-01 2.746e-02 4.703 0.000240 ***
## BOOK -4.336e-03 -1.885e-02 1.215e-02 -0.357 0.725883
## log_ADD_FEE 2.568e+00 7.431e-02 2.010e+00 1.277 0.219712
## R_B_COST -3.258e-03 -1.042e-01 2.189e-03 -1.488 0.156128
## SF_RATIO -1.300e+00 -1.876e-01 5.055e-01 -2.572 0.020479 *
## GRADUAT 1.590e-01 8.799e-02 1.657e-01 0.960 0.351347
## OUT_STAT 1.837e-03 1.893e-01 1.079e-03 1.704 0.107803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.753 on 16 degrees of freedom
## (29 пропущенных наблюдений удалены)
## Multiple R-squared: 0.9616, Adjusted R-squared: 0.9448
## F-statistic: 57.28 on 7 and 16 DF, p-value: 3.83e-10
AIC(model_pr_noPhD)
## [1] 168.0601
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)
## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.5973852 0.30419786
## BOOK BOOK 0.3840774 -0.08885404
## R_B_COST R_B_COST 0.7551732 -0.34870910
## OUT_STAT OUT_STAT 0.9141040 0.39183509
## AVRCOMB AVRCOMB 0.9692789 0.76171345
## SF_RATIO SF_RATIO 0.8252002 -0.54081295
## GRADUAT GRADUAT 0.8544490 0.23336928
Попробуем убрать BOOK (из-за низкой частной корреляции) и GRADUAT
model_pr1<-lm(NEW10 ~ AVRCOMB + log_ADD_FEE + R_B_COST + SF_RATIO + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr1)
##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + log_ADD_FEE + R_B_COST + SF_RATIO +
## OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.4607 -1.9482 -0.7267 2.3280 13.4017
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.185e+02 NA 2.745e+01 -4.315 0.000417 ***
## AVRCOMB 1.372e-01 6.446e-01 2.323e-02 5.907 1.36e-05 ***
## log_ADD_FEE 2.961e+00 8.567e-02 1.880e+00 1.575 0.132718
## R_B_COST -3.886e-03 -1.243e-01 2.040e-03 -1.905 0.072923 .
## SF_RATIO -1.336e+00 -1.927e-01 4.914e-01 -2.718 0.014095 *
## OUT_STAT 2.164e-03 2.229e-01 9.870e-04 2.192 0.041756 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.583 on 18 degrees of freedom
## (29 пропущенных наблюдений удалены)
## Multiple R-squared: 0.959, Adjusted R-squared: 0.9476
## F-statistic: 84.15 on 5 and 18 DF, p-value: 7.676e-12
AIC(model_pr1)
## [1] 165.6631
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D", "BOOK", "GRADUAT"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)
## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.5684235 0.3479827
## R_B_COST R_B_COST 0.7447506 -0.4095647
## OUT_STAT OUT_STAT 0.9088596 0.4590401
## AVRCOMB AVRCOMB 0.9668878 0.8122173
## SF_RATIO SF_RATIO 0.8236772 -0.5394659
AIC(model_pr_new)
## [1] 148.3376
AIC(model_pr_noPhD)
## [1] 168.0601
AIC(model_pr1)
## [1] 165.6631
library(olsrr)
##
## Присоединяю пакет: 'olsrr'
## Следующий объект скрыт от 'package:MASS':
##
## cement
## Следующий объект скрыт от 'package:datasets':
##
## rivers
backward_aic<-ols_step_backward_aic(model_pr_new, progress=TRUE, details=TRUE)
## Backward Elimination Method
## ---------------------------
##
## Candidate Terms:
##
## 1. AVRCOMB
## 2. BOOK
## 3. log_ADD_FEE
## 4. R_B_COST
## 5. PH_D
## 6. SF_RATIO
## 7. GRADUAT
## 8. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 148.3376
##
## Initiating stepwise selection...
##
## Table: Removing Existing Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## BOOK 1 146.734 156.135 97.586 0.96397 0.94457
## PH_D 1 147.994 157.395 97.513 0.96174 0.94114
## GRADUAT 1 148.658 158.058 97.491 0.96051 0.93925
## R_B_COST 1 149.412 158.813 97.480 0.95907 0.93703
## log_ADD_FEE 1 149.693 159.094 97.480 0.95852 0.93618
## OUT_STAT 1 150.834 160.235 97.501 0.95620 0.93262
## SF_RATIO 1 153.877 163.277 97.737 0.94937 0.92211
## AVRCOMB 1 156.737 166.137 98.197 0.94199 0.91075
## -----------------------------------------------------------------------
##
## Step => 1
## Removed => BOOK
## Model => NEW10 ~ AVRCOMB + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 146.7343
##
## Table: Removing Existing Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## PH_D 1 147.170 155.526 94.423 0.95954 0.94220
## GRADUAT 1 147.861 156.217 94.521 0.95818 0.94026
## log_ADD_FEE 1 147.962 156.319 94.536 0.95798 0.93997
## R_B_COST 1 148.265 156.622 94.583 0.95737 0.93910
## OUT_STAT 1 150.933 159.289 95.079 0.95160 0.93085
## SF_RATIO 1 153.017 161.373 95.578 0.94655 0.92364
## AVRCOMB 1 155.056 163.412 96.162 0.94110 0.91585
## -----------------------------------------------------------------------
##
##
## No more variables to be removed.
##
## Variables Removed:
##
## => BOOK
forward_aic<-ols_step_forward_aic(model_pr_new, progress=TRUE, details=TRUE)
## Forward Selection Method
## ------------------------
##
## Candidate Terms:
##
## 1. AVRCOMB
## 2. BOOK
## 3. log_ADD_FEE
## 4. R_B_COST
## 5. PH_D
## 6. SF_RATIO
## 7. GRADUAT
## 8. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ 1
## AIC => 202.5251
##
## Initiating stepwise selection...
##
## Table: Adding New Variables
## -------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## AVRCOMB 1 147.617 150.751 87.731 0.93346 0.92995
## OUT_STAT 1 171.142 174.276 107.789 0.79601 0.78527
## GRADUAT 1 176.803 179.937 112.954 0.73289 0.71883
## PH_D 1 186.429 189.563 121.963 0.57756 0.55533
## SF_RATIO 1 189.574 192.707 124.955 0.50933 0.48350
## log_ADD_FEE 1 202.476 205.609 137.417 0.09297 0.04523
## R_B_COST 1 203.410 206.543 138.328 0.05173 0.00182
## BOOK 1 204.456 207.590 139.350 0.00327 -0.04919
## -------------------------------------------------------------------------
##
## Step => 1
## Added => AVRCOMB
## Model => NEW10 ~ AVRCOMB
## AIC => 147.617
##
## Table: Adding New Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## SF_RATIO 1 145.049 149.227 86.340 0.94647 0.94052
## OUT_STAT 1 147.493 151.671 88.069 0.93986 0.93318
## BOOK 1 147.975 152.153 88.412 0.93846 0.93162
## log_ADD_FEE 1 148.581 152.759 88.846 0.93666 0.92962
## PH_D 1 149.095 153.274 89.215 0.93509 0.92788
## GRADUAT 1 149.166 153.344 89.265 0.93487 0.92764
## R_B_COST 1 149.275 153.453 89.344 0.93453 0.92726
## -----------------------------------------------------------------------
##
## Step => 2
## Added => SF_RATIO
## Model => NEW10 ~ AVRCOMB + SF_RATIO
## AIC => 145.0489
##
## Table: Adding New Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## BOOK 1 145.266 150.488 87.603 0.95082 0.94215
## OUT_STAT 1 145.397 150.620 87.679 0.95052 0.94178
## GRADUAT 1 146.151 151.373 88.121 0.94871 0.93966
## log_ADD_FEE 1 146.479 151.701 88.314 0.94790 0.93871
## R_B_COST 1 146.891 152.113 88.557 0.94687 0.93749
## PH_D 1 146.968 152.190 88.603 0.94667 0.93726
## -----------------------------------------------------------------------
##
##
## No more variables to be added.
##
## Variables Selected:
##
## => AVRCOMB
## => SF_RATIO
Получили, что forward_aic имеет лучшее значение, чем backward_aic, то есть, является более хорошей линейной моделью. К тому же, у forward_aic меньше признаков.
forward_aic
##
##
## Stepwise Summary
## -------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## 0 Base Model 202.525 204.614 139.293 0.00000 0.00000
## 1 AVRCOMB 147.617 150.751 87.731 0.93346 0.92995
## 2 SF_RATIO 145.049 149.227 86.340 0.94647 0.94052
## -------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ---------------------------------------------------------------
## R 0.973 RMSE 6.323
## R-Squared 0.946 MSE 39.975
## Adj. R-Squared 0.941 Coef. Var 11.813
## Pred R-Squared 0.925 AIC 145.049
## MAE 4.914 SBC 149.227
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ----------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ----------------------------------------------------------------------
## Regression 14841.759 2 7420.879 159.117 0.0000
## Residual 839.479 18 46.638
## Total 15681.238 20
## ----------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------------
## (Intercept) -150.963 21.538 -7.009 0.000 -196.212 -105.714
## AVRCOMB 0.188 0.015 0.869 12.124 0.000 0.155 0.220
## SF_RATIO -0.986 0.471 -0.150 -2.091 0.051 -1.977 0.004
## ----------------------------------------------------------------------------------------------
backward_aic
##
##
## Stepwise Summary
## -------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## 0 Full Model 148.338 158.783 101.117 0.96464 0.94107
## 1 BOOK 146.734 156.135 96.228 0.96397 0.94457
## -------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ---------------------------------------------------------------
## R 0.982 RMSE 5.187
## R-Squared 0.964 MSE 26.905
## Adj. R-Squared 0.945 Coef. Var 11.404
## Pred R-Squared 0.901 AIC 146.734
## MAE 4.081 SBC 156.135
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ---------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ---------------------------------------------------------------------
## Regression 15116.223 7 2159.460 49.685 0.0000
## Residual 565.015 13 43.463
## Total 15681.238 20
## ---------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------------
## (Intercept) -71.623 48.490 -1.477 0.163 -176.379 33.133
## AVRCOMB 0.110 0.038 0.511 2.873 0.013 0.027 0.193
## log_ADD_FEE 3.162 2.151 0.097 1.470 0.165 -1.486 7.810
## R_B_COST -0.004 0.002 -0.126 -1.543 0.147 -0.009 0.001
## PH_D -0.573 0.453 -0.147 -1.264 0.228 -1.552 0.406
## SF_RATIO -1.413 0.564 -0.215 -2.507 0.026 -2.631 -0.195
## GRADUAT 0.344 0.238 0.201 1.445 0.172 -0.171 0.860
## OUT_STAT 0.003 0.001 0.286 2.113 0.055 0.000 0.005
## -------------------------------------------------------------------------------------------
Здесь можно сравнить по R^2, например.
model_pr_new<-lm(NEW10 ~ AVRCOMB + SF_RATIO, data=pr_data)|>lm.beta()
summary(model_pr_new)
##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + SF_RATIO, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.0779 -6.4604 0.1359 4.5401 17.2921
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -133.57926 NA 20.45055 -6.532 7.66e-07 ***
## AVRCOMB 0.17283 0.84064 0.01434 12.055 6.52e-12 ***
## SF_RATIO -1.18802 -0.17520 0.47288 -2.512 0.0188 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.591 on 25 degrees of freedom
## (25 пропущенных наблюдений удалены)
## Multiple R-squared: 0.9309, Adjusted R-squared: 0.9254
## F-statistic: 168.4 on 2 and 25 DF, p-value: 3.104e-15
AIC(model_pr_new)
## [1] 197.7974
shapiro.test(residuals(model_pr_new))
##
## Shapiro-Wilk normality test
##
## data: residuals(model_pr_new)
## W = 0.9525, p-value = 0.2288
plot(model_pr_new, 1:5)
У 15-го сильное плечо (ещё видно, как линия идёт к ней)
plot(residuals(model_pr_new), rstudent(model_pr_new))
Все значения не так сильно отличаются, так что посмотрим на outliers по Куку и Махаланобису
# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pr_new)
cook_df <- data.frame(
obs = seq_along(cooks_vals),
cook = cooks_vals
) %>%
arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))
plot(
x, cook_df$cook,
type = "h",
main = "Cook's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Cook's distance",
xaxt = "n" # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)
# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pr_new)[,-1],
center = colMeans(model.matrix(model_pr_new)[,-1]),
cov = cov(model.matrix(model_pr_new)[,-1]))
mah_df <- data.frame(
obs = seq_along(mahalanobis_distance),
mah = mahalanobis_distance
) %>%
arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05
# Scatter plot для расстояния Махаланобиса
plot(
x, mah_df$mah,
type = "h",
main = "Mahalanobis's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Mahalanobis's distance",
xaxt = "n"
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)
Не нашёл outliers (в случае AVRCOMB)
set.seed(456) # Другое зерно для генерации новых значений
max_AVRCOMB<-max(pr_data$AVRCOMB, na.rm=TRUE)
max_SF_RATIO<-max(pr_data$SF_RATIO, na.rm=TRUE)
new_data <- data.frame(
AVRCOMB = max_AVRCOMB + 1:8,
SF_RATIO = max_SF_RATIO + 1:8
)
# Прогнозируем значение NEW10
predicted_values <- predict(model_pr_new, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")
## Прогнозируемые значения NEW10:
print(predicted_values)
## 1 2 3 4 5 6 7 8
## 84.73871 83.72351 82.70832 81.69312 80.67793 79.66273 78.64754 77.63234
# 6. Доверительные и предсказательные интервалы
# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pr_new, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")
## [1] "Доверительные интервалы (95%):\n"
print(confidence_intervals)
## fit lwr upr
## 1 84.73871 68.04653 101.4309
## 2 83.72351 66.08992 101.3571
## 3 82.70832 64.12789 101.2887
## 4 81.69312 62.16121 101.2250
## 5 80.67793 60.19053 101.1653
## 6 79.66273 58.21640 101.1091
## 7 78.64754 56.23924 101.0558
## 8 77.63234 54.25944 101.0052
# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pr_new, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")
## [1] "Предсказательные интервалы (95%):\n"
print(prediction_intervals)
## fit lwr upr
## 1 84.73871 61.86841 107.6090
## 2 83.72351 60.15733 107.2897
## 3 82.70832 58.42552 106.9911
## 4 81.69312 56.67478 106.7115
## 5 80.67793 54.90672 106.4491
## 6 79.66273 53.12282 106.2026
## 7 78.64754 51.32441 105.9707
## 8 77.63234 49.51270 105.7520
# Создаем data.frame для прогнозов и интервалов
predictions_df <- data.frame(
AVRCOMB = new_data$AVRCOMB,
SF_RATIO = new_data$SF_RATIO,
Predicted = predicted_values,
Lower_Conf = confidence_intervals[, "lwr"],
Upper_Conf = confidence_intervals[, "upr"],
Lower_Pred = prediction_intervals[, "lwr"],
Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df
# Визуализация (только для демонстрации, нужен более сложный график для нескольких предикторов)
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) + # Условная ось x, т.к. несколько предикторов
geom_point() +
geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
x = "Прогноз (номер)",
y = "NEW10") +
scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
theme_bw()
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.9, label = "Синие: Доверительные интервалы", color = "blue", hjust = 0) +
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.8, label = "Красные: Предсказательные интервалы", color = "red", hjust = 0)
pub_data<-col_I_regr_split$public
print(paste0("Число наблюдений: ", dim(pub_data)[1]))
## [1] "Число наблюдений: 123"
colSums(is.na(pub_data))
## ...1 log_ADD_FEE BOOK R_B_COST OUT_STAT NEW10
## 0 33 2 0 1 15
## AVRCOMB SF_RATIO PH_D GRADUAT PPIND
## 43 0 5 4 0
Много пропусков в AVRCOMB
plot_numeric_pairs(pub_data, "Общественные вузы")
outlier_pub_out_eye<-pub_data[pub_data$OUT_STAT>14001,]
outlier_pub_out_eye
plot_numeric_pairs(pub_data, "Общественные вузы", highlight = list(red=outlier_pub_out_eye))
Здесь условия из .tsk файла
Стоит помнить, что у AVRCOMB много NA. В такой модели не будет выбросов “на глаз”
model_pub_new <- lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT,
data = pub_data)|>lm.beta()
summary(model_pub_new)
##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST +
## PH_D + SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.462 -7.954 -1.811 5.800 33.454
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.957e+02 NA 5.205e+01 -3.760 0.000531 ***
## AVRCOMB 1.124e-01 3.828e-01 4.176e-02 2.691 0.010254 *
## BOOK 3.210e-02 1.521e-01 2.409e-02 1.333 0.190033
## log_ADD_FEE 4.645e+00 2.355e-01 2.416e+00 1.922 0.061536 .
## R_B_COST 4.859e-03 1.906e-01 3.400e-03 1.429 0.160534
## PH_D 4.937e-01 1.328e-01 4.894e-01 1.009 0.318999
## SF_RATIO 7.921e-02 1.507e-02 6.456e-01 0.123 0.902959
## GRADUAT 1.657e-01 1.096e-01 2.419e-01 0.685 0.497348
## OUT_STAT -8.391e-04 -9.592e-02 1.337e-03 -0.628 0.533709
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.48 on 41 degrees of freedom
## (73 пропущенных наблюдений удалены)
## Multiple R-squared: 0.5774, Adjusted R-squared: 0.4949
## F-statistic: 7.002 on 8 and 41 DF, p-value: 8.678e-06
plot_confidence_ellipse(model_pub_new, c(3, 2), level=0.95)
На примере этой пары: AVRCOMB и SF_RATIO значимы
regressors <- setdiff(names(pub_data), c("...1", "PPIND", "NEW10"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pub_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pub_data[,regressors_full])
)
print(redundancy_table)
## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.6083961 0.2875337
## BOOK BOOK 0.4915332 0.2037467
## R_B_COST R_B_COST 0.6691166 0.2178369
## OUT_STAT OUT_STAT 0.7502290 -0.0975559
## AVRCOMB AVRCOMB 0.7530387 0.3874765
## SF_RATIO SF_RATIO 0.5634165 0.0191561
## PH_D PH_D 0.6475179 0.1556248
## GRADUAT GRADUAT 0.7761031 0.1063354
Можем убрать OUT_STAT и GRADUAT
AIC(model_pub_new)
## [1] 425.9309
model_pub1 <- lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO,
data = pub_data)|>lm.beta()
summary(model_pub1)
##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST +
## PH_D + SF_RATIO, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.863 -7.446 -2.303 7.075 32.472
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -2.035e+02 NA 4.107e+01 -4.955 1.07e-05 ***
## AVRCOMB 1.151e-01 4.081e-01 3.088e-02 3.727 0.000539 ***
## BOOK 3.398e-02 1.619e-01 2.134e-02 1.592 0.118345
## log_ADD_FEE 5.201e+00 2.632e-01 2.298e+00 2.263 0.028499 *
## R_B_COST 3.886e-03 1.546e-01 2.844e-03 1.367 0.178554
## PH_D 5.734e-01 1.550e-01 4.594e-01 1.248 0.218431
## SF_RATIO 1.030e-01 1.985e-02 5.417e-01 0.190 0.849989
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.14 on 45 degrees of freedom
## (71 пропущенное наблюдение удалено)
## Multiple R-squared: 0.5601, Adjusted R-squared: 0.5014
## F-statistic: 9.548 on 6 and 45 DF, p-value: 9.133e-07
AIC(model_pub1)
## [1] 438.6636
AIC увеличился, а значит предпочтительнее выбирать model_pub_new.
regr_high <- setdiff(names(pub_data), c("...1", "PPIND", "NEW10", "OUT_STAT", "GRADUAT"))
regr_fullhigh<-c(regr_high, "NEW10")
redund_high <- data.frame(
Регрессор = regr_high,
Множественная_корреляция = sapply(regr_high, calc_multicollinearity, pub_data[, regr_fullhigh]),
Частная_корреляция = sapply(regr_high, calc_partial_correlation, pub_data[,regr_fullhigh])
)
print(redund_high)
## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.5924855 0.31966785
## BOOK BOOK 0.3246989 0.23093274
## R_B_COST R_B_COST 0.5163902 0.19961448
## AVRCOMB AVRCOMB 0.6139954 0.48569728
## SF_RATIO SF_RATIO 0.3206426 0.02834558
## PH_D PH_D 0.6222610 0.18292398
library(olsrr)
backward_aic_pub<-ols_step_backward_aic(model_pub_new, progress=TRUE, details=TRUE)
## Backward Elimination Method
## ---------------------------
##
## Candidate Terms:
##
## 1. AVRCOMB
## 2. BOOK
## 3. log_ADD_FEE
## 4. R_B_COST
## 5. PH_D
## 6. SF_RATIO
## 7. GRADUAT
## 8. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 425.9309
##
## Initiating stepwise selection...
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## SF_RATIO 1 423.949 441.157 285.464 0.57725 0.50679
## OUT_STAT 1 424.409 441.617 285.755 0.57334 0.50223
## GRADUAT 1 424.499 441.708 285.813 0.57257 0.50133
## PH_D 1 425.157 442.365 286.230 0.56691 0.49473
## BOOK 1 426.051 443.259 286.802 0.55910 0.48562
## R_B_COST 1 426.362 443.570 287.002 0.55635 0.48241
## log_ADD_FEE 1 428.246 445.454 288.223 0.53932 0.46253
## AVRCOMB 1 432.065 449.273 290.751 0.50275 0.41987
## ------------------------------------------------------------------------
##
## Step => 1
## Removed => SF_RATIO
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + GRADUAT + OUT_STAT
## AIC => 423.9492
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## GRADUAT 1 422.552 437.849 283.437 0.57212 0.51241
## OUT_STAT 1 422.682 437.978 283.525 0.57101 0.51115
## PH_D 1 423.231 438.527 283.901 0.56627 0.50575
## BOOK 1 424.053 439.349 284.466 0.55908 0.49755
## R_B_COST 1 424.548 439.844 284.807 0.55470 0.49256
## log_ADD_FEE 1 426.379 441.676 286.079 0.53808 0.47363
## AVRCOMB 1 430.065 445.361 288.678 0.50275 0.43336
## ------------------------------------------------------------------------
##
## Step => 2
## Removed => GRADUAT
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + OUT_STAT
## AIC => 422.5524
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## OUT_STAT 1 420.913 434.297 281.292 0.56902 0.52004
## PH_D 1 421.997 435.381 282.087 0.55957 0.50952
## BOOK 1 422.446 435.830 282.417 0.55560 0.50510
## R_B_COST 1 423.127 436.511 282.920 0.54951 0.49831
## log_ADD_FEE 1 425.924 439.308 284.998 0.52358 0.46945
## AVRCOMB 1 435.531 448.915 292.326 0.42266 0.35706
## ------------------------------------------------------------------------
##
## Step => 3
## Removed => OUT_STAT
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D
## AIC => 420.9132
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## PH_D 1 420.513 431.985 280.150 0.55501 0.51545
## R_B_COST 1 421.160 432.632 280.657 0.54921 0.50914
## BOOK 1 421.831 433.303 281.185 0.54312 0.50250
## log_ADD_FEE 1 424.596 436.068 283.366 0.51715 0.47423
## AVRCOMB 1 433.721 445.193 290.696 0.42046 0.36895
## ------------------------------------------------------------------------
##
## Step => 4
## Removed => PH_D
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST
## AIC => 420.5128
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## R_B_COST 1 421.792 431.352 280.561 0.52484 0.49386
## BOOK 1 422.106 431.666 280.822 0.52185 0.49067
## log_ADD_FEE 1 425.733 435.293 283.858 0.48588 0.45235
## AVRCOMB 1 436.611 446.171 293.093 0.36093 0.31925
## ------------------------------------------------------------------------
##
##
## No more variables to be removed.
##
## Variables Removed:
##
## => SF_RATIO
## => GRADUAT
## => OUT_STAT
## => PH_D
forward_aic_pub<-ols_step_forward_aic(model_pub_new, progress=TRUE, details=TRUE)
## Forward Selection Method
## ------------------------
##
## Candidate Terms:
##
## 1. AVRCOMB
## 2. BOOK
## 3. log_ADD_FEE
## 4. R_B_COST
## 5. PH_D
## 6. SF_RATIO
## 7. GRADUAT
## 8. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ 1
## AIC => 452.9976
##
## Initiating stepwise selection...
##
## Table: Adding New Variables
## -------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## AVRCOMB 1 432.770 438.506 290.015 0.35889 0.34553
## PH_D 1 440.325 446.061 297.005 0.25432 0.23878
## log_ADD_FEE 1 440.669 446.405 297.324 0.24917 0.23353
## GRADUAT 1 440.775 446.511 297.423 0.24757 0.23189
## R_B_COST 1 445.964 451.700 302.247 0.16529 0.14790
## BOOK 1 450.887 456.623 306.843 0.07893 0.05974
## OUT_STAT 1 454.461 460.197 310.192 0.01067 -0.00994
## SF_RATIO 1 454.634 460.370 310.354 0.00724 -0.01344
## -------------------------------------------------------------------------
##
## Step => 1
## Added => AVRCOMB
## Model => NEW10 ~ AVRCOMB
## AIC => 432.7703
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## log_ADD_FEE 1 423.754 431.402 281.872 0.48566 0.46377
## R_B_COST 1 426.389 434.038 284.194 0.45782 0.43475
## PH_D 1 427.487 435.135 285.163 0.44579 0.42221
## BOOK 1 431.895 439.543 289.066 0.39471 0.36895
## GRADUAT 1 433.164 440.812 290.193 0.37916 0.35274
## OUT_STAT 1 434.743 442.391 291.599 0.35923 0.33197
## SF_RATIO 1 434.747 442.395 291.602 0.35919 0.33192
## ------------------------------------------------------------------------
##
## Step => 2
## Added => log_ADD_FEE
## Model => NEW10 ~ AVRCOMB + log_ADD_FEE
## AIC => 423.7541
##
## Table: Adding New Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## BOOK 1 421.792 431.352 280.561 0.52484 0.49386
## R_B_COST 1 422.106 431.666 280.822 0.52185 0.49067
## PH_D 1 422.181 431.741 280.885 0.52113 0.48990
## OUT_STAT 1 425.569 435.129 283.721 0.48756 0.45414
## GRADUAT 1 425.579 435.140 283.729 0.48745 0.45403
## SF_RATIO 1 425.746 435.306 283.869 0.48574 0.45221
## ----------------------------------------------------------------------
##
## Step => 3
## Added => BOOK
## Model => NEW10 ~ AVRCOMB + log_ADD_FEE + BOOK
## AIC => 421.7919
##
## Table: Adding New Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## R_B_COST 1 420.513 431.985 280.150 0.55501 0.51545
## PH_D 1 421.160 432.632 280.657 0.54921 0.50914
## GRADUAT 1 422.984 434.456 282.092 0.53246 0.49091
## OUT_STAT 1 423.748 435.220 282.695 0.52526 0.48306
## SF_RATIO 1 423.783 435.255 282.723 0.52493 0.48270
## ----------------------------------------------------------------------
##
## Step => 4
## Added => R_B_COST
## Model => NEW10 ~ AVRCOMB + log_ADD_FEE + BOOK + R_B_COST
## AIC => 420.5128
##
## Table: Adding New Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## PH_D 1 420.913 434.297 281.292 0.56902 0.52004
## OUT_STAT 1 421.997 435.381 282.087 0.55957 0.50952
## GRADUAT 1 422.228 435.613 282.257 0.55753 0.50725
## SF_RATIO 1 422.498 435.882 282.455 0.55514 0.50459
## ----------------------------------------------------------------------
##
##
## No more variables to be added.
##
## Variables Selected:
##
## => AVRCOMB
## => log_ADD_FEE
## => BOOK
## => R_B_COST
Получили совпадение результатов пошагового forward и пошагового backward.
forward_aic_pub
##
##
## Stepwise Summary
## --------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## --------------------------------------------------------------------------
## 0 Base Model 452.998 456.822 309.665 0.00000 0.00000
## 1 AVRCOMB 432.770 438.506 290.015 0.35889 0.34553
## 2 log_ADD_FEE 423.754 431.402 281.872 0.48566 0.46377
## 3 BOOK 421.792 431.352 280.561 0.52484 0.49386
## 4 R_B_COST 420.513 431.985 280.150 0.55501 0.51545
## --------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ----------------------------------------------------------------
## R 0.745 RMSE 14.385
## R-Squared 0.555 MSE 206.930
## Adj. R-Squared 0.515 Coef. Var 42.498
## Pred R-Squared 0.424 AIC 420.513
## MAE 11.265 SBC 431.985
## ----------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ---------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ---------------------------------------------------------------------
## Regression 12904.403 4 3226.101 14.031 0.0000
## Residual 10346.477 45 229.922
## Total 23250.880 49
## ---------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------------
## (Intercept) -184.284 31.816 -5.792 0.000 -248.364 -120.203
## AVRCOMB 0.136 0.031 0.462 4.430 0.000 0.074 0.197
## log_ADD_FEE 5.845 2.211 0.296 2.644 0.011 1.392 10.297
## BOOK 0.039 0.021 0.185 1.831 0.074 -0.004 0.082
## R_B_COST 0.005 0.003 0.192 1.746 0.088 -0.001 0.011
## ----------------------------------------------------------------------------------------------
backward_aic_pub
##
##
## Stepwise Summary
## -------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## 0 Full Model 425.931 445.051 287.892 0.57740 0.49494
## 1 SF_RATIO 423.949 441.157 285.030 0.57725 0.50679
## 2 GRADUAT 422.552 437.849 282.885 0.57212 0.51241
## 3 OUT_STAT 420.913 434.297 280.618 0.56902 0.52004
## 4 PH_D 420.513 431.985 279.705 0.55501 0.51545
## -------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ----------------------------------------------------------------
## R 0.745 RMSE 14.385
## R-Squared 0.555 MSE 206.930
## Adj. R-Squared 0.515 Coef. Var 42.498
## Pred R-Squared 0.424 AIC 420.513
## MAE 11.265 SBC 431.985
## ----------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ---------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ---------------------------------------------------------------------
## Regression 12904.403 4 3226.101 14.031 0.0000
## Residual 10346.477 45 229.922
## Total 23250.880 49
## ---------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------------
## (Intercept) -184.284 31.816 -5.792 0.000 -248.364 -120.203
## AVRCOMB 0.136 0.031 0.462 4.430 0.000 0.074 0.197
## BOOK 0.039 0.021 0.185 1.831 0.074 -0.004 0.082
## log_ADD_FEE 5.845 2.211 0.296 2.644 0.011 1.392 10.297
## R_B_COST 0.005 0.003 0.192 1.746 0.088 -0.001 0.011
## ----------------------------------------------------------------------------------------------
model_pub_new<-lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST, data=pub_data)|>lm.beta()
summary(model_pub_new)
##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST,
## data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.246 -9.548 -1.513 9.069 29.958
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.895e+02 NA 2.957e+01 -6.410 5.02e-08 ***
## AVRCOMB 1.261e-01 4.023e-01 2.927e-02 4.308 7.70e-05 ***
## BOOK 4.740e-02 2.089e-01 2.025e-02 2.341 0.0233 *
## log_ADD_FEE 6.998e+00 3.365e-01 2.135e+00 3.278 0.0019 **
## R_B_COST 5.708e-03 2.171e-01 2.666e-03 2.141 0.0371 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.23 on 50 degrees of freedom
## (68 пропущенных наблюдений удалены)
## Multiple R-squared: 0.6167, Adjusted R-squared: 0.5861
## F-statistic: 20.11 on 4 and 50 DF, p-value: 6.357e-10
AIC(model_pub_new)
## [1] 462.4016
shapiro.test(residuals(model_pub_new))
##
## Shapiro-Wilk normality test
##
## data: residuals(model_pub_new)
## W = 0.97633, p-value = 0.3475
plot(model_pub_new)
На первом графике Residuals vs. Fitted заметен сильный изгиб (матожидание ошибки ненулевое)
plot(residuals(model_pub_new), rstudent(model_pub_new))
Все значения не так сильно отличаются, так что посмотрим на outliers по Куку и Махаланобису
# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pub_new)
cook_df <- data.frame(
obs = seq_along(cooks_vals),
cook = cooks_vals
) %>%
arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))
plot(
x, cook_df$cook,
type = "h",
main = "Cook's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Cook's distance",
xaxt = "n" # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)
# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pub_new)[,-1],
center = colMeans(model.matrix(model_pub_new)[,-1]),
cov = cov(model.matrix(model_pub_new)[,-1]))
mah_df <- data.frame(
obs = seq_along(mahalanobis_distance),
mah = mahalanobis_distance
) %>%
arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05
# Scatter plot для расстояния Махаланобиса
plot(
x, mah_df$mah,
type = "h",
main = "Mahalanobis's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Mahalanobis's distance",
xaxt = "n"
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)
Не нашёл outliers (в случае AVRCOMB)
set.seed(456) # Другое зерно для генерации новых значений
max_AVRCOMB<-max(pub_data$AVRCOMB, na.rm=TRUE)
max_BOOK<-max(pub_data$BOOK, na.rm=TRUE)
max_log_ADD_FEE<-max(pub_data$log_ADD_FEE, na.rm=TRUE)
max_R_B_COST<-max(pub_data$R_B_COST, na.rm=TRUE)
new_data <- data.frame(
AVRCOMB = max_AVRCOMB - seq(50, 120, by = 10),
BOOK = max_BOOK - seq(10, 80, by = 10),
log_ADD_FEE = max_log_ADD_FEE - seq(0.01, 0.08, by=0.01),
R_B_COST = max_R_B_COST - seq(10, 80, by=10)
)
# Прогнозируем значение NEW10
predicted_values <- predict(model_pub_new, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")
## Прогнозируемые значения NEW10:
print(predicted_values)
## 1 2 3 4 5 6 7 8
## 96.97503 95.11299 93.25095 91.38891 89.52687 87.66483 85.80280 83.94076
# 6. Доверительные и предсказательные интервалы
# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pub_new, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")
## [1] "Доверительные интервалы (95%):\n"
print(confidence_intervals)
## fit lwr upr
## 1 96.97503 81.95568 111.99437
## 2 95.11299 80.48533 109.74065
## 3 93.25095 78.99581 107.50609
## 4 91.38891 77.48558 105.29224
## 5 89.52687 75.95304 103.10071
## 6 87.66483 74.39651 100.93315
## 7 85.80280 72.81432 98.79127
## 8 83.94076 71.20476 96.67675
# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pub_new, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")
## [1] "Предсказательные интервалы (95%):\n"
print(prediction_intervals)
## fit lwr upr
## 1 96.97503 62.89617 131.0539
## 2 95.11299 61.20493 129.0210
## 3 93.25095 59.50192 127.0000
## 4 91.38891 57.78697 124.9909
## 5 89.52687 56.05992 122.9938
## 6 87.66483 54.32063 121.0090
## 7 85.80280 52.56895 119.0366
## 8 83.94076 50.80477 117.0767
predictions_df <- data.frame(
AVRCOMB = new_data$AVRCOMB,
BOOK = new_data$BOOK,
log_ADD_FEE = new_data$log_ADD_FEE,
R_B_COST = new_data$R_B_COST,
Predicted = predicted_values,
Lower_Conf = confidence_intervals[, "lwr"],
Upper_Conf = confidence_intervals[, "upr"],
Lower_Pred = prediction_intervals[, "lwr"],
Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) +
geom_point() +
geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
x = "Прогноз (номер)",
y = "NEW10") +
scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
theme_bw()
Важные пункты для исправления:
Использование AVRCOMB (особенно в случае частных вузов) не интересно. Поэтому уберём её из рассматриваемых признаков.
Доверительный эллипсоид некорректно строится, поэтому нужно использовать другую функцию построения доверительного эллипсоида.
Рассмотреть прологарифмированный NEW10 для общественных вузов
Рассматриваемые наблюдения: 1. PPIND (для группировки) 2. ADD_FEE - дополнительные выплаты 3. BOOK - плата за бронирование 4. NEW10 (зависимая переменная, также NEW25) - студентов из топ 10% своей старшей школы 5. AVRCOMB - средний комбинированный балл (лишний признак) 6. SF_RATIO - Student/Faculty ratio 7. PH_D - Pct. of faculty with Ph.D.’s 8. GRADUAT - Graduation rate 9. R_B_COST - Room and board costs 10. OUT_STAT - Out-of-state tuition
colNames<-c("...1", "ADD_FEE", "BOOK", "R_B_COST","OUT_STAT","NEW10", "SF_RATIO","PH_D","GRADUAT","PPIND")
col_I_regr<-col_I_sn[c(colNames)]
col_I_regr<-col_I_regr|>mutate(ADD_FEE = log(ADD_FEE))|>rename(log_ADD_FEE=ADD_FEE)
col_I_regr_split<-split(col_I_regr, col_I_regr$PPIND)
# Автор: -Николай-
library(GGally)
library(ggplot2)
library(dplyr)
library(rlang)
plot_numeric_pairs <- function(df, title = NULL, highlight = NULL) {
# Выбираем имена всех numeric-столбцов
numeric_cols <- df %>% dplyr::select(dplyr::where(is.numeric)) %>% names()
# Функция для нижних панелей с возможностью выделения точек разными цветами
custom_points <- function(data, mapping, ...) {
# Отрисовываем базовые точки синим цветом
p <- ggplot(data, mapping) +
geom_point(alpha = 0.6, size = 2, color = "blue", ...)
if (!is.null(highlight)) {
# Если highlight не является списком, превращаем его в список с одним элементом и цветом "red"
if (!is.list(highlight) || inherits(highlight, "data.frame")) {
highlight <- list(red = highlight)
}
# Получаем имена переменных из mapping с помощью as_label()
xvar <- as_label(mapping$x)
yvar <- as_label(mapping$y)
# Для каждого набора точек в списке highlight
for (col in names(highlight)) {
hl_df <- highlight[[col]]
# Проверяем, что hl_df содержит нужные переменные
if(all(c(xvar, yvar) %in% names(hl_df))) {
# Находим совпадающие строки между данными панели и hl_df
data_highlight <- semi_join(data, hl_df, by = c(xvar, yvar))
# Накладываем выделенные точки нужного цвета
p <- p + geom_point(data = data_highlight, mapping = mapping,
alpha = 0.6, size = 2, color = col, ...)
}
}
}
p
}
# Строим матрицу парных графиков
ggpairs(df,
columns = numeric_cols,
columnLabels = numeric_cols,
title = title,
upper = list(continuous = wrap("cor", size = 3)),
lower = list(continuous = custom_points),
diag = list(continuous = wrap("barDiag", bins = 10, fill = "lightblue", color = "black"))
) +
theme_bw() +
theme(
strip.text.x = element_text(size = 7, angle = 0, hjust = 0),
strip.text.y = element_text(size = 8, angle = 0, hjust = 0),
axis.text = element_text(size = 6),
plot.margin = unit(c(1, 2, 1, 1), "cm")
)
}
ggpairs(col_I_regr[,-which(names(col_I_regr)=="...1")], aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Комментарий: замечаем, что для публичных вузов NEW10 имеет выраженный правый хвост. Поэтому стоит использовать log_NEW10 для публичных вузов.
При постановке вопроса о построении линейной модели с категориальным признаком стоит обратить внимание на то, близки ли корреляции разных групп. Если они отличаются, то лучше рассматривать каждую категорию по отдельности. Если близки, то вместе. Это связано с тем, что разные корреляции по разному влияют на общую модель, будут возникать неоднородности.
pr_data<-col_I_regr_split$private
print(paste0("Число наблюдений: ", dim(pr_data)[1]))
## [1] "Число наблюдений: 53"
colSums(is.na(pr_data))
## ...1 log_ADD_FEE BOOK R_B_COST OUT_STAT NEW10
## 0 7 0 0 0 1
## SF_RATIO PH_D GRADUAT PPIND
## 0 4 1 0
plot_numeric_pairs(pr_data, "Частные вузы")
outlier_pr_book_eye<-pr_data[pr_data$BOOK>875,]
outlier_pr_out_eye<-pr_data[pr_data$OUT_STAT<5000,]
outlier_pr_book_eye
outlier_pr_out_eye
plot_numeric_pairs(pr_data, "Частные вузы", list(red = outlier_pr_book_eye, green = outlier_pr_out_eye))
library(lm.beta)
model_pr_full <- lm(NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT,
data = pr_data) |> lm.beta()
summary(model_pr_full)
##
## Call:
## lm(formula = NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO +
## GRADUAT + OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.1636 -6.8538 0.0723 7.2414 25.5012
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -9.289705 NA 38.856630 -0.239 0.81252
## BOOK 0.002300 0.015211 0.013944 0.165 0.86999
## log_ADD_FEE 1.877765 0.063416 2.978846 0.630 0.53280
## R_B_COST -0.002897 -0.120242 0.002971 -0.975 0.33658
## PH_D -0.219550 -0.069638 0.508017 -0.432 0.66843
## SF_RATIO -1.772400 -0.324548 0.604579 -2.932 0.00608 **
## GRADUAT 1.069124 0.642349 0.215538 4.960 2.07e-05 ***
## OUT_STAT 0.001308 0.192473 0.001098 1.192 0.24180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.4 on 33 degrees of freedom
## (12 пропущенных наблюдений удалены)
## Multiple R-squared: 0.73, Adjusted R-squared: 0.6727
## F-statistic: 12.75 on 7 and 33 DF, p-value: 8.508e-08
plot_confidence_ellipse(model_pr_full, c(2, 3))
plot_confidence_ellipse(model_pr_full, c(4, 5))
plot_confidence_ellipse(model_pr_full, c(6, 7))
plot_confidence_ellipse(model_pr_full, c(8, 2))
# Автор: -Евгений-
# Вычисление множественной корреляции
calc_multicollinearity <- function(var, data) {
X <- as.matrix(data[, setdiff(names(data), var)]) # Матрица остальных регрессоров
y <- as.matrix(data[[var]]) # Текущий регрессор как зависимая переменная
model <- lm(y ~ X)
sqrt(summary(model)$r.squared)
}
calc_partial_correlation <- function(var, data) {
# Остатки регрессии текущего регрессора на остальные регрессоры
data<-data|>na.omit()
model_X <- lm(data[[var]] ~ ., data = data[, !names(data) %in% c(var, "NEW10")])
res_X <- residuals(model_X)
# Остатки регрессии Y на остальные регрессоры
model_Y <- lm(NEW10 ~ ., data = data[, !names(data) %in% var])
res_Y <- residuals(model_Y)
cor(res_X, res_Y, use = "complete.obs") # Корреляция остатков
}
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)
## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.4485260 0.10907796
## BOOK BOOK 0.1965671 0.02870252
## R_B_COST R_B_COST 0.6906122 -0.16735732
## OUT_STAT OUT_STAT 0.8361740 0.20315778
## SF_RATIO SF_RATIO 0.6858132 -0.45456011
## PH_D PH_D 0.8286465 -0.07501937
## GRADUAT GRADUAT 0.8488216 0.65354784
Комментарий: Стоит исключить PH_D из-за высокой множественной и низкой частной корреляций (также можно исключить BOOK из-за низкой частной корреляции).
AIC(model_pr_full)
## [1] 338.2841
summary(model_pr_full)
##
## Call:
## lm(formula = NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO +
## GRADUAT + OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.1636 -6.8538 0.0723 7.2414 25.5012
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -9.289705 NA 38.856630 -0.239 0.81252
## BOOK 0.002300 0.015211 0.013944 0.165 0.86999
## log_ADD_FEE 1.877765 0.063416 2.978846 0.630 0.53280
## R_B_COST -0.002897 -0.120242 0.002971 -0.975 0.33658
## PH_D -0.219550 -0.069638 0.508017 -0.432 0.66843
## SF_RATIO -1.772400 -0.324548 0.604579 -2.932 0.00608 **
## GRADUAT 1.069124 0.642349 0.215538 4.960 2.07e-05 ***
## OUT_STAT 0.001308 0.192473 0.001098 1.192 0.24180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.4 on 33 degrees of freedom
## (12 пропущенных наблюдений удалены)
## Multiple R-squared: 0.73, Adjusted R-squared: 0.6727
## F-statistic: 12.75 on 7 and 33 DF, p-value: 8.508e-08
model_pr_nono<-lm(NEW10 ~ log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr_nono)
##
## Call:
## lm(formula = NEW10 ~ log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT +
## OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.1896 -7.8512 0.4492 6.5775 29.4475
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.254e+01 NA 2.304e+01 -0.544 0.5895
## log_ADD_FEE 1.313e+00 4.221e-02 2.810e+00 0.467 0.6430
## R_B_COST -3.726e-03 -1.462e-01 2.743e-03 -1.358 0.1823
## SF_RATIO -1.872e+00 -3.275e-01 5.701e-01 -3.284 0.0022 **
## GRADUAT 1.019e+00 5.877e-01 1.880e-01 5.420 3.56e-06 ***
## OUT_STAT 1.228e-03 1.710e-01 9.368e-04 1.311 0.1976
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.49 on 38 degrees of freedom
## (9 пропущенных наблюдений удалены)
## Multiple R-squared: 0.726, Adjusted R-squared: 0.6899
## F-statistic: 20.13 on 5 and 38 DF, p-value: 9.083e-10
AIC(model_pr_nono)
## [1] 361.386
Комментарий: замечаем, что после того, как убрали BOOK и PH_D увеличился Adj.R^2, немного уменьшился R^2, уменьшился p-value, но увеличился AIC.
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D", "BOOK"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)
## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.3480743 0.07558012
## R_B_COST R_B_COST 0.6376041 -0.21521296
## OUT_STAT OUT_STAT 0.7708944 0.20809285
## SF_RATIO SF_RATIO 0.6597111 -0.47020518
## GRADUAT GRADUAT 0.8087224 0.66027976
Попробуем убрать log_ADD_FEE
model_pr1<-lm(NEW10 ~ GRADUAT + R_B_COST + SF_RATIO + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr1)
##
## Call:
## lm(formula = NEW10 ~ GRADUAT + R_B_COST + SF_RATIO + OUT_STAT,
## data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.041 -9.947 -1.776 8.355 49.615
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 14.1924930 NA 21.7789414 0.652 0.51786
## GRADUAT 0.9698293 0.6097925 0.1999565 4.850 1.45e-05 ***
## R_B_COST -0.0038847 -0.1556376 0.0029142 -1.333 0.18908
## SF_RATIO -1.6934700 -0.3080606 0.6201252 -2.731 0.00893 **
## OUT_STAT 0.0003467 0.0543364 0.0010089 0.344 0.73268
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.36 on 46 degrees of freedom
## (2 пропущенных наблюдений удалены)
## Multiple R-squared: 0.6291, Adjusted R-squared: 0.5968
## F-statistic: 19.5 on 4 and 46 DF, p-value: 1.918e-09
AIC(model_pr1)
## [1] 430.1013
Комментарий: AIC увеличилась, а R^2 сильно уменьшились, p-value увеличился.
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D", "BOOK", "log_ADD_FEE"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)
## Регрессор Множественная_корреляция Частная_корреляция
## R_B_COST R_B_COST 0.6561089 -0.19285565
## OUT_STAT OUT_STAT 0.8235937 0.05060282
## SF_RATIO SF_RATIO 0.6743579 -0.37350259
## GRADUAT GRADUAT 0.8139325 0.58168914
AIC(model_pr_full)
## [1] 338.2841
AIC(model_pr_nono)
## [1] 361.386
AIC(model_pr1)
## [1] 430.1013
library(olsrr)
backward_aic<-ols_step_backward_aic(model_pr_full, progress=TRUE, details=TRUE)
## Backward Elimination Method
## ---------------------------
##
## Candidate Terms:
##
## 1. BOOK
## 2. log_ADD_FEE
## 3. R_B_COST
## 4. PH_D
## 5. SF_RATIO
## 6. GRADUAT
## 7. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 338.2841
##
## Initiating stepwise selection...
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## BOOK 1 336.318 350.026 223.228 0.72977 0.68209
## PH_D 1 336.516 350.224 223.348 0.72847 0.68055
## log_ADD_FEE 1 336.775 350.483 223.505 0.72674 0.67852
## R_B_COST 1 337.449 351.157 223.916 0.72221 0.67319
## OUT_STAT 1 338.012 351.721 224.262 0.71837 0.66867
## SF_RATIO 1 345.774 359.483 229.220 0.65968 0.59962
## GRADUAT 1 359.125 372.833 238.570 0.52868 0.44551
## ------------------------------------------------------------------------
##
## Step => 1
## Removed => BOOK
## Model => NEW10 ~ log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 336.3179
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## PH_D 1 334.535 346.530 220.890 0.72834 0.68953
## log_ADD_FEE 1 334.821 346.816 221.080 0.72644 0.68736
## R_B_COST 1 335.500 347.495 221.536 0.72187 0.68213
## OUT_STAT 1 336.045 348.040 221.902 0.71815 0.67789
## SF_RATIO 1 343.786 355.781 227.258 0.65958 0.61094
## GRADUAT 1 357.284 369.279 237.270 0.52684 0.45925
## ------------------------------------------------------------------------
##
## Step => 2
## Removed => PH_D
## Model => NEW10 ~ log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 334.5351
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## log_ADD_FEE 1 332.891 343.172 218.679 0.72597 0.69552
## R_B_COST 1 333.513 343.794 219.132 0.72178 0.69087
## OUT_STAT 1 334.111 344.392 219.570 0.71770 0.68633
## SF_RATIO 1 341.927 352.208 225.394 0.65840 0.62045
## GRADUAT 1 358.480 368.761 238.448 0.48849 0.43166
## ------------------------------------------------------------------------
##
## Step => 3
## Removed => log_ADD_FEE
## Model => NEW10 ~ R_B_COST + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 332.891
##
## Table: Removing Existing Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## R_B_COST 1 331.621 340.189 216.806 0.72105 0.69843
## OUT_STAT 1 332.295 340.863 217.339 0.71642 0.69343
## SF_RATIO 1 339.927 348.495 223.430 0.65840 0.63071
## GRADUAT 1 358.030 366.598 238.458 0.46878 0.42571
## ----------------------------------------------------------------------
##
## Step => 4
## Removed => R_B_COST
## Model => NEW10 ~ SF_RATIO + GRADUAT + OUT_STAT
## AIC => 331.6207
##
## Table: Removing Existing Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## OUT_STAT 1 330.381 337.236 215.046 0.71582 0.70087
## SF_RATIO 1 338.121 344.975 221.632 0.65678 0.63872
## GRADUAT 1 358.218 365.073 239.135 0.43966 0.41017
## ----------------------------------------------------------------------
##
## Step => 5
## Removed => OUT_STAT
## Model => NEW10 ~ SF_RATIO + GRADUAT
## AIC => 330.3813
##
## Table: Removing Existing Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## SF_RATIO 1 339.980 345.121 223.161 0.62291 0.61324
## GRADUAT 1 365.570 370.711 246.575 0.29611 0.27806
## ----------------------------------------------------------------------
##
##
## No more variables to be removed.
##
## Variables Removed:
##
## => BOOK
## => PH_D
## => log_ADD_FEE
## => R_B_COST
## => OUT_STAT
backward_aic
##
##
## Stepwise Summary
## --------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## --------------------------------------------------------------------------
## 0 Full Model 338.284 353.706 225.692 0.73000 0.67272
## 1 BOOK 336.318 350.026 222.763 0.72977 0.68209
## 2 PH_D 334.535 346.530 220.181 0.72834 0.68953
## 3 log_ADD_FEE 332.891 343.172 217.888 0.72597 0.69552
## 4 R_B_COST 331.621 340.189 216.109 0.72105 0.69843
## 5 OUT_STAT 330.381 337.236 214.490 0.71582 0.70087
## --------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ----------------------------------------------------------------
## R 0.846 RMSE 12.336
## R-Squared 0.716 MSE 152.183
## Adj. R-Squared 0.701 Coef. Var 23.175
## Pred R-Squared 0.672 AIC 330.381
## MAE 10.318 SBC 337.236
## ----------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ---------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ---------------------------------------------------------------------
## Regression 15717.001 2 7858.501 47.86 0.0000
## Residual 6239.487 38 164.197
## Total 21956.488 40
## ---------------------------------------------------------------------
##
## Parameter Estimates
## ------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ------------------------------------------------------------------------------------------
## (Intercept) -18.639 14.637 -1.273 0.211 -48.270 10.992
## SF_RATIO -1.760 0.499 -0.322 -3.525 0.001 -2.770 -0.749
## GRADUAT 1.140 0.152 0.685 7.492 0.000 0.832 1.448
## ------------------------------------------------------------------------------------------
forward_aic<-ols_step_forward_aic(model_pr_full, progress=TRUE, details=TRUE)
## Forward Selection Method
## ------------------------
##
## Candidate Terms:
##
## 1. BOOK
## 2. log_ADD_FEE
## 3. R_B_COST
## 4. PH_D
## 5. SF_RATIO
## 6. GRADUAT
## 7. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ 1
## AIC => 377.966
##
## Initiating stepwise selection...
##
## Table: Adding New Variables
## -------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## GRADUAT 1 339.980 345.121 223.161 0.62291 0.61324
## PH_D 1 358.045 363.186 239.617 0.41413 0.39910
## OUT_STAT 1 362.074 367.215 243.334 0.35363 0.33706
## SF_RATIO 1 365.570 370.711 246.575 0.29611 0.27806
## R_B_COST 1 378.009 383.150 258.223 0.04661 0.02216
## BOOK 1 379.594 384.735 259.720 0.00903 -0.01638
## log_ADD_FEE 1 379.787 384.928 259.903 0.00435 -0.02118
## -------------------------------------------------------------------------
##
## Step => 1
## Added => GRADUAT
## Model => NEW10 ~ GRADUAT
## AIC => 339.9804
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## SF_RATIO 1 330.381 337.236 215.046 0.71582 0.70087
## OUT_STAT 1 338.121 344.975 221.632 0.65678 0.63872
## PH_D 1 340.187 347.041 223.400 0.63905 0.62005
## R_B_COST 1 341.456 348.310 224.490 0.62770 0.60810
## log_ADD_FEE 1 341.939 348.793 224.905 0.62329 0.60346
## BOOK 1 341.973 348.827 224.934 0.62297 0.60313
## ------------------------------------------------------------------------
##
## Step => 2
## Added => SF_RATIO
## Model => NEW10 ~ GRADUAT + SF_RATIO
## AIC => 330.3813
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## OUT_STAT 1 331.621 340.189 216.806 0.72105 0.69843
## log_ADD_FEE 1 332.265 340.833 217.315 0.71663 0.69366
## PH_D 1 332.272 340.840 217.320 0.71658 0.69360
## R_B_COST 1 332.295 340.863 217.339 0.71642 0.69343
## BOOK 1 332.321 340.889 217.359 0.71624 0.69323
## ------------------------------------------------------------------------
##
##
## No more variables to be added.
##
## Variables Selected:
##
## => GRADUAT
## => SF_RATIO
forward_aic
##
##
## Stepwise Summary
## -------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## 0 Base Model 377.966 381.393 259.401 0.00000 0.00000
## 1 GRADUAT 339.980 345.121 223.161 0.62291 0.61324
## 2 SF_RATIO 330.381 337.236 215.046 0.71582 0.70087
## -------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ----------------------------------------------------------------
## R 0.846 RMSE 12.336
## R-Squared 0.716 MSE 152.183
## Adj. R-Squared 0.701 Coef. Var 23.175
## Pred R-Squared 0.672 AIC 330.381
## MAE 10.318 SBC 337.236
## ----------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ---------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ---------------------------------------------------------------------
## Regression 15717.001 2 7858.501 47.86 0.0000
## Residual 6239.487 38 164.197
## Total 21956.488 40
## ---------------------------------------------------------------------
##
## Parameter Estimates
## ------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ------------------------------------------------------------------------------------------
## (Intercept) -18.639 14.637 -1.273 0.211 -48.270 10.992
## GRADUAT 1.140 0.152 0.685 7.492 0.000 0.832 1.448
## SF_RATIO -1.760 0.499 -0.322 -3.525 0.001 -2.770 -0.749
## ------------------------------------------------------------------------------------------
Комментарий: строим Backward и Forward AIC. Backward AIC исключил всё, кроме SF_RATIO и GRADUAT. Forward AIC включил только GRADUAT и SF_RATIO. Оба пошаговых алгоритма дали один и тот же результат.
Интерпретация параметров: число лучших 10% студентов NEW10 увеличивается с увеличением выпускников GRADUAT и уменьшается с числом студентов на преподавателя SF_RATIO.
Здесь можно сравнить по R^2, например.
model_pr_full<-lm(NEW10 ~ GRADUAT + SF_RATIO, data=pr_data)
summary(model_pr_full)
##
## Call:
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.763 -8.954 -2.873 9.404 52.149
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.1676 16.4411 -0.314 0.75465
## GRADUAT 0.9740 0.1626 5.990 2.6e-07 ***
## SF_RATIO -1.5182 0.5621 -2.701 0.00952 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.35 on 48 degrees of freedom
## (2 пропущенных наблюдений удалены)
## Multiple R-squared: 0.6133, Adjusted R-squared: 0.5972
## F-statistic: 38.07 on 2 and 48 DF, p-value: 1.247e-10
plot(model_pr_full, 1:5)
У 51-го сильное плечо (ещё видно, как линия идёт к ней)
plot(residuals(model_pr_full), rstudent(model_pr_full))
Комментарий: замечаем точку (51), которая сильно влияет на линейную модель.
influence.measures(model_pr_full)
## Influence measures of
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data) :
##
## dfb.1_ dfb.GRAD dfb.SF_R dffit cov.r cook.d hat inf
## 1 0.156252 -0.045588 -2.35e-01 0.4007 0.816 4.95e-02 0.0317
## 2 -0.142405 0.171118 7.76e-02 0.2199 1.057 1.61e-02 0.0498
## 3 0.007595 -0.007708 -9.15e-05 0.0145 1.099 7.16e-05 0.0308
## 4 0.007442 0.011854 -7.56e-02 -0.1148 1.115 4.46e-03 0.0579
## 5 -0.062742 0.118154 -4.21e-02 0.2030 1.079 1.38e-02 0.0554
## 6 0.209642 -0.134914 -3.80e-01 -0.4260 1.077 5.98e-02 0.1036
## 7 -0.010436 0.018846 -3.62e-02 -0.0813 1.093 2.24e-03 0.0360
## 8 -0.156587 0.114409 8.70e-02 -0.3109 0.849 3.03e-02 0.0231
## 9 -0.157795 0.122821 1.25e-01 -0.2056 1.037 1.41e-02 0.0387
## 10 0.025289 -0.043555 6.54e-03 -0.0741 1.102 1.86e-03 0.0411
## 12 -0.031867 0.027892 2.43e-02 -0.0335 1.206 3.81e-04 0.1170 *
## 13 0.007036 -0.029448 3.02e-02 -0.0739 1.117 1.85e-03 0.0523
## 14 0.148220 -0.160722 -2.28e-02 0.2043 1.101 1.40e-02 0.0676
## 15 -0.440485 0.369605 3.53e-01 -0.4841 0.980 7.54e-02 0.0802
## 16 -0.020974 0.058786 -3.49e-02 0.1383 1.069 6.44e-03 0.0357
## 17 -0.013268 -0.014281 4.90e-02 -0.0892 1.099 2.70e-03 0.0421
## 18 -0.165762 0.166156 1.42e-01 0.2044 1.139 1.41e-02 0.0899
## 19 -0.037008 0.028235 2.21e-02 -0.0639 1.082 1.39e-03 0.0251
## 20 -0.014143 0.000378 3.18e-02 -0.0433 1.136 6.37e-04 0.0640
## 21 -0.169351 0.147972 1.93e-01 0.2218 1.183 1.66e-02 0.1200
## 22 0.024521 -0.044146 -1.97e-02 -0.1356 1.034 6.15e-03 0.0219
## 24 -0.233928 0.276537 1.44e-01 0.3683 0.934 4.35e-02 0.0451
## 25 -0.081062 0.088996 1.72e-02 -0.1024 1.203 3.56e-03 0.1191 *
## 26 0.071188 -0.084065 -4.77e-02 -0.1185 1.086 4.74e-03 0.0399
## 27 -0.182509 0.145727 1.67e-01 -0.2044 1.152 1.41e-02 0.0979
## 28 -0.067432 -0.008354 1.65e-01 -0.2386 1.063 1.90e-02 0.0563
## 29 -0.008973 0.028263 -2.49e-02 0.0639 1.127 1.39e-03 0.0585
## 30 -0.110949 0.144068 3.66e-02 0.1853 1.086 1.15e-02 0.0551
## 31 0.084879 -0.059446 -1.54e-01 -0.1872 1.099 1.18e-02 0.0624
## 32 0.022119 -0.040726 1.33e-02 -0.0689 1.121 1.61e-03 0.0550
## 33 -0.016792 0.037689 -1.59e-02 0.0784 1.095 2.09e-03 0.0372
## 34 0.154239 -0.125010 -2.45e-01 -0.3179 0.989 3.30e-02 0.0484
## 35 -0.052304 0.065696 -2.44e-02 -0.1143 1.110 4.43e-03 0.0545
## 36 0.189408 -0.151498 -1.43e-01 0.2417 1.014 1.93e-02 0.0395
## 37 -0.029891 -0.001702 1.42e-01 0.2076 1.068 1.44e-02 0.0514
## 38 -0.105057 0.106146 9.90e-03 -0.1824 1.035 1.11e-02 0.0329
## 39 -0.060391 0.032813 7.43e-02 -0.1045 1.091 3.70e-03 0.0397
## 40 -0.018011 0.057624 -5.00e-02 0.1320 1.105 5.89e-03 0.0543
## 41 0.373278 -0.291574 -3.66e-01 0.4214 1.143 5.90e-02 0.1335
## 42 0.041672 -0.031904 -2.39e-02 0.0735 1.077 1.83e-03 0.0247
## 43 -0.091131 0.106561 -1.98e-02 -0.1722 1.081 9.97e-03 0.0494
## 44 0.288062 -0.292022 -2.68e-01 -0.4058 0.948 5.29e-02 0.0556
## 45 -0.015808 0.050662 -3.76e-02 0.1224 1.083 5.06e-03 0.0394
## 46 -0.062436 0.092305 4.33e-03 0.1352 1.091 6.17e-03 0.0466
## 47 0.024276 -0.009007 -3.60e-02 0.0538 1.102 9.82e-04 0.0377
## 48 -0.073400 0.113328 1.46e-02 0.1985 1.019 1.31e-02 0.0315
## 49 0.000713 0.003235 -2.22e-02 -0.0423 1.099 6.09e-04 0.0335
## 50 -0.004980 0.008760 -1.19e-02 -0.0260 1.121 2.29e-04 0.0506
## 51 1.078464 -1.600184 8.12e-01 2.5174 0.450 1.48e+00 0.2292 *
## 52 0.047770 -0.034990 -1.09e-01 -0.1676 1.051 9.41e-03 0.0348
## 53 0.027573 -0.015210 -5.76e-02 -0.0660 1.184 1.48e-03 0.1022
row.names(model_pr_full$model)
## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "12" "13" "14" "15" "16"
## [16] "17" "18" "19" "20" "21" "22" "24" "25" "26" "27" "28" "29" "30" "31" "32"
## [31] "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45" "46" "47"
## [46] "48" "49" "50" "51" "52" "53"
outlier_cook<-pr_data|>dplyr::filter(NEW10==48&SF_RATIO==20.5)
outlier_cook
Комментарий: выбросом оказался Brigham Young University.
# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pr_full)
cook_df <- data.frame(
obs = row.names(model_pr_full$model),
cook = cooks_vals
) %>%
arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))
plot(
x, cook_df$cook,
type = "h",
main = "Cook's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Cook's distance",
xaxt = "n" # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)
Комментарий: высокие значения у 29 и 38.
# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pr_full)[,-1],
center = colMeans(model.matrix(model_pr_full)[,-1]),
cov = cov(model.matrix(model_pr_full)[,-1]))
mah_df <- data.frame(
obs = row.names(model_pr_full$model),
mah = mahalanobis_distance
) %>%
arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05
# Scatter plot для расстояния Махаланобиса
plot(
x, mah_df$mah,
type = "h",
main = "Mahalanobis's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Mahalanobis's distance",
xaxt = "n"
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)
Комментарий: 29-ое наблюдение имеет большое расстояние по Махаланобису.
pr_data_outlier<-pr_data[-c(51),]
plot_numeric_pairs(pr_data_outlier)
model_pr_outlier<-lm(NEW10~GRADUAT + SF_RATIO, data=pr_data_outlier)
summary(model_pr_full|>lm.beta())
##
## Call:
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.763 -8.954 -2.873 9.404 52.149
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -5.1676 NA 16.4411 -0.314 0.75465
## GRADUAT 0.9740 0.6124 0.1626 5.990 2.6e-07 ***
## SF_RATIO -1.5182 -0.2762 0.5621 -2.701 0.00952 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.35 on 48 degrees of freedom
## (2 пропущенных наблюдений удалены)
## Multiple R-squared: 0.6133, Adjusted R-squared: 0.5972
## F-statistic: 38.07 on 2 and 48 DF, p-value: 1.247e-10
summary(model_pr_outlier|>lm.beta())
##
## Call:
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data_outlier)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.6995 -8.8741 0.3934 8.0222 31.4411
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -20.0312 NA 14.1534 -1.415 0.163576
## GRADUAT 1.1921 0.6758 0.1443 8.263 1.03e-10 ***
## SF_RATIO -1.9009 -0.3250 0.4784 -3.973 0.000242 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.87 on 47 degrees of freedom
## (2 пропущенных наблюдений удалены)
## Multiple R-squared: 0.7332, Adjusted R-squared: 0.7218
## F-statistic: 64.57 on 2 and 47 DF, p-value: 3.284e-14
Комментарий: замечаем существенную разницу в R^2, p-value и значимости SF_RATIO.
plot(model_pr_outlier, 1:5)
AIC(model_pr_full)
## [1] 428.221
AIC(model_pr_outlier)
## [1] 402.2787
Можем заметить, что модель без выброса лучше по AIC, чем с выбросом.
set.seed(456) # Другое зерно для генерации новых значений
max_GRADUAT<-max(pr_data$GRADUAT, na.rm=TRUE)
max_SF_RATIO<-max(pr_data$SF_RATIO, na.rm=TRUE)
new_data <- data.frame(
GRADUAT = max_GRADUAT - 1:8,
SF_RATIO = max_SF_RATIO - 1:8
)
# Прогнозируем значение NEW10
predicted_values <- predict(model_pr_outlier, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")
## Прогнозируемые значения NEW10:
print(predicted_values)
## 1 2 3 4 5 6 7 8
## 59.72933 60.43809 61.14684 61.85560 62.56435 63.27311 63.98186 64.69061
# 6. Доверительные и предсказательные интервалы
# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pr_outlier, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")
## [1] "Доверительные интервалы (95%):\n"
print(confidence_intervals)
## fit lwr upr
## 1 59.72933 46.77615 72.68252
## 2 60.43809 48.53109 72.34508
## 3 61.14684 50.27377 72.01991
## 4 61.85560 52.00033 71.71086
## 5 62.56435 53.70520 71.42350
## 6 63.27311 55.38017 71.16604
## 7 63.98186 57.01279 70.95093
## 8 64.69061 58.58381 70.79742
# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pr_outlier, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")
## [1] "Предсказательные интервалы (95%):\n"
print(prediction_intervals)
## fit lwr upr
## 1 59.72933 30.78140 88.67727
## 2 60.43809 31.94292 88.93325
## 3 61.14684 33.06800 89.22569
## 4 61.85560 34.15498 89.55621
## 5 62.56435 35.20229 89.92641
## 6 63.27311 36.20843 90.33778
## 7 63.98186 37.17205 90.79167
## 8 64.69061 38.09191 91.28932
# Создаем data.frame для прогнозов и интервалов
predictions_df <- data.frame(
GRADUAT = new_data$GRADUAT,
SF_RATIO = new_data$SF_RATIO,
Predicted = predicted_values,
Lower_Conf = confidence_intervals[, "lwr"],
Upper_Conf = confidence_intervals[, "upr"],
Lower_Pred = prediction_intervals[, "lwr"],
Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df
# Визуализация (только для демонстрации, нужен более сложный график для нескольких предикторов)
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) + # Условная ось x, т.к. несколько предикторов
geom_point() +
geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
x = "Прогноз (номер)",
y = "NEW10") +
scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
theme_bw()
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.9, label = "Синие: Доверительные интервалы", color = "blue", hjust = 0) +
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.8, label = "Красные: Предсказательные интервалы", color = "red", hjust = 0)
pub_data<-col_I_regr_split$public
pub_data<-pub_data|>mutate(NEW10 = log(NEW10))|>rename(log_NEW10 = NEW10)
print(paste0("Число наблюдений: ", dim(pub_data)[1]))
## [1] "Число наблюдений: 123"
colSums(is.na(pub_data))
## ...1 log_ADD_FEE BOOK R_B_COST OUT_STAT log_NEW10
## 0 33 2 0 1 15
## SF_RATIO PH_D GRADUAT PPIND
## 0 5 4 0
Много пропусков по log_ADD_FEE (около четверти).
plot_numeric_pairs(pub_data, "Общественные вузы")
outlier_pub_out_eye<-pub_data[pub_data$OUT_STAT>14000,]
outlier_pub_out_eye
plot_numeric_pairs(pub_data, "Общественные вузы", list(red = outlier_pub_out_eye))
library(lm.beta)
model_pub_full <- lm(log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT,
data = pub_data)
summary(model_pub_full|>lm.beta())
##
## Call:
## lm(formula = log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D +
## SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23717 -0.22303 0.04705 0.22529 0.89431
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -4.821e-02 NA 8.271e-01 -0.058 0.9537
## BOOK 8.739e-04 1.632e-01 5.489e-04 1.592 0.1161
## log_ADD_FEE 4.078e-03 7.717e-03 5.368e-02 0.076 0.9397
## R_B_COST 8.749e-05 1.289e-01 7.275e-05 1.203 0.2334
## PH_D 2.303e-02 2.706e-01 8.974e-03 2.566 0.0126 *
## SF_RATIO -9.484e-03 -7.112e-02 1.289e-02 -0.736 0.4645
## GRADUAT 1.970e-02 5.464e-01 4.132e-03 4.768 1.07e-05 ***
## OUT_STAT -5.986e-05 -2.567e-01 2.867e-05 -2.088 0.0407 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4196 on 66 degrees of freedom
## (49 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4814, Adjusted R-squared: 0.4264
## F-statistic: 8.752 on 7 and 66 DF, p-value: 1.449e-07
plot_confidence_ellipse(model_pub_full, c(2, 3))
plot_confidence_ellipse(model_pub_full, c(4, 5))
plot_confidence_ellipse(model_pub_full, c(6, 7))
plot_confidence_ellipse(model_pub_full, c(8, 2))
# Автор: -Евгений-
# Вычисление множественной корреляции
calc_multicollinearity <- function(var, data) {
X <- as.matrix(data[, setdiff(names(data), var)]) # Матрица остальных регрессоров
y <- as.matrix(data[[var]]) # Текущий регрессор как зависимая переменная
model <- lm(y ~ X)
sqrt(summary(model)$r.squared)
}
calc_partial_correlation <- function(var, data) {
# Остатки регрессии текущего регрессора на остальные регрессоры
data<-data|>na.omit()
model_X <- lm(data[[var]] ~ ., data = data[, !names(data) %in% c(var, "log_NEW10")])
res_X <- residuals(model_X)
# Остатки регрессии Y на остальные регрессоры
model_Y <- lm(log_NEW10 ~ ., data = data[, !names(data) %in% var])
res_Y <- residuals(model_Y)
cor(res_X, res_Y, use = "complete.obs") # Корреляция остатков
}
regressors <- setdiff(names(pub_data), c("...1", "PPIND", "log_NEW10"))
regressors_full<-c(regressors, "log_NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pub_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pub_data[,regressors_full])
)
print(redundancy_table)
## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.4883582 0.009351778
## BOOK BOOK 0.5291391 0.192313348
## R_B_COST R_B_COST 0.5747925 0.146434350
## OUT_STAT OUT_STAT 0.7158033 -0.248923976
## SF_RATIO SF_RATIO 0.4073614 -0.090198497
## PH_D PH_D 0.5981262 0.301176912
## GRADUAT GRADUAT 0.7449868 0.506130402
Исключим log_ADD_FEE и SF_RATIO (так как у них низкие частные корреляции)
AIC(model_pub_full)
## [1] 90.99499
summary(model_pub_full)
##
## Call:
## lm(formula = log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D +
## SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23717 -0.22303 0.04705 0.22529 0.89431
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.821e-02 8.271e-01 -0.058 0.9537
## BOOK 8.739e-04 5.489e-04 1.592 0.1161
## log_ADD_FEE 4.078e-03 5.368e-02 0.076 0.9397
## R_B_COST 8.749e-05 7.275e-05 1.203 0.2334
## PH_D 2.303e-02 8.974e-03 2.566 0.0126 *
## SF_RATIO -9.484e-03 1.289e-02 -0.736 0.4645
## GRADUAT 1.970e-02 4.132e-03 4.768 1.07e-05 ***
## OUT_STAT -5.986e-05 2.867e-05 -2.088 0.0407 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4196 on 66 degrees of freedom
## (49 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4814, Adjusted R-squared: 0.4264
## F-statistic: 8.752 on 7 and 66 DF, p-value: 1.449e-07
model_pub_nono<-lm(log_NEW10 ~ R_B_COST + PH_D + GRADUAT + OUT_STAT + BOOK, data=pub_data)
summary(model_pub_full|>lm.beta())
##
## Call:
## lm(formula = log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D +
## SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23717 -0.22303 0.04705 0.22529 0.89431
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -4.821e-02 NA 8.271e-01 -0.058 0.9537
## BOOK 8.739e-04 1.632e-01 5.489e-04 1.592 0.1161
## log_ADD_FEE 4.078e-03 7.717e-03 5.368e-02 0.076 0.9397
## R_B_COST 8.749e-05 1.289e-01 7.275e-05 1.203 0.2334
## PH_D 2.303e-02 2.706e-01 8.974e-03 2.566 0.0126 *
## SF_RATIO -9.484e-03 -7.112e-02 1.289e-02 -0.736 0.4645
## GRADUAT 1.970e-02 5.464e-01 4.132e-03 4.768 1.07e-05 ***
## OUT_STAT -5.986e-05 -2.567e-01 2.867e-05 -2.088 0.0407 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4196 on 66 degrees of freedom
## (49 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4814, Adjusted R-squared: 0.4264
## F-statistic: 8.752 on 7 and 66 DF, p-value: 1.449e-07
AIC(model_pub_full)
## [1] 90.99499
summary(model_pub_nono|>lm.beta())
##
## Call:
## lm(formula = log_NEW10 ~ R_B_COST + PH_D + GRADUAT + OUT_STAT +
## BOOK, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.21996 -0.24955 -0.00251 0.24842 0.81900
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 7.390e-02 NA 5.501e-01 0.134 0.8934
## R_B_COST 1.272e-04 1.806e-01 6.475e-05 1.964 0.0525 .
## PH_D 1.748e-02 2.156e-01 7.462e-03 2.342 0.0213 *
## GRADUAT 1.808e-02 4.905e-01 3.432e-03 5.268 9.03e-07 ***
## OUT_STAT -5.027e-05 -2.003e-01 2.389e-05 -2.104 0.0381 *
## BOOK 1.044e-03 1.903e-01 4.510e-04 2.315 0.0228 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4225 on 92 degrees of freedom
## (25 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4772, Adjusted R-squared: 0.4488
## F-statistic: 16.8 on 5 and 92 DF, p-value: 9.211e-12
AIC(model_pub_nono)
## [1] 117.0556
Комментарий: AIC увеличился, p-value уменьшился. Теперь каждый признак, кроме R_B_COST является значимым (при alpha=0.05).
regressors <- setdiff(names(pub_data), c("...1", "PPIND", "log_NEW10", "log_ADD_FEE", "SF_RATIO"))
regressors_full<-c(regressors, "log_NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pub_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pub_data[,regressors_full])
)
print(redundancy_table)
## Регрессор Множественная_корреляция Частная_корреляция
## BOOK BOOK 0.4535409 0.2346185
## R_B_COST R_B_COST 0.5958627 0.2006126
## OUT_STAT OUT_STAT 0.6338016 -0.2142731
## PH_D PH_D 0.6059169 0.2372086
## GRADUAT GRADUAT 0.7046461 0.4813910
Тут уже не так понятно. Перейдём к пошаговому AIC.
AIC(model_pub_full)
## [1] 90.99499
AIC(model_pub_nono)
## [1] 117.0556
library(olsrr)
backward_aic<-ols_step_backward_aic(model_pub_full, progress=TRUE, details=TRUE)
## Backward Elimination Method
## ---------------------------
##
## Candidate Terms:
##
## 1. BOOK
## 2. log_ADD_FEE
## 3. R_B_COST
## 4. PH_D
## 5. SF_RATIO
## 6. GRADUAT
## 7. OUT_STAT
##
##
## Step => 0
## Model => log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 90.99499
##
## Initiating stepwise selection...
##
## Table: Removing Existing Variables
## -------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## log_ADD_FEE 1 89.001 107.434 -119.335 0.48135 0.43490
## SF_RATIO 1 89.599 108.032 -118.859 0.47714 0.43032
## R_B_COST 1 90.599 109.032 -118.062 0.47003 0.42257
## BOOK 1 91.784 110.216 -117.115 0.46148 0.41325
## OUT_STAT 1 93.728 112.161 -115.555 0.44714 0.39763
## PH_D 1 96.032 114.464 -113.699 0.42966 0.37858
## GRADUAT 1 110.895 129.327 -101.488 0.30279 0.24035
## -------------------------------------------------------------------------
##
## Step => 1
## Removed => log_ADD_FEE
## Model => log_NEW10 ~ BOOK + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 89.00146
##
## Table: Removing Existing Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## SF_RATIO 1 87.604 103.733 -121.080 0.47711 0.43866
## R_B_COST 1 88.775 104.904 -120.111 0.46877 0.42970
## BOOK 1 89.835 105.964 -119.232 0.46110 0.42148
## OUT_STAT 1 91.997 108.125 -117.436 0.44513 0.40433
## PH_D 1 94.185 110.314 -115.609 0.42847 0.38645
## GRADUAT 1 111.437 127.565 -100.972 0.27842 0.22537
## -----------------------------------------------------------------------
##
## Step => 2
## Removed => SF_RATIO
## Model => log_NEW10 ~ BOOK + R_B_COST + PH_D + GRADUAT + OUT_STAT
## AIC => 87.60406
##
## Table: Removing Existing Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## R_B_COST 1 87.268 101.092 -121.878 0.46522 0.43422
## BOOK 1 88.652 102.476 -120.689 0.45512 0.42354
## OUT_STAT 1 90.007 103.831 -119.523 0.44505 0.41288
## PH_D 1 94.291 108.115 -115.824 0.41198 0.37789
## GRADUAT 1 109.486 123.311 -102.540 0.27794 0.23608
## -----------------------------------------------------------------------
##
## Step => 3
## Removed => R_B_COST
## Model => log_NEW10 ~ BOOK + PH_D + GRADUAT + OUT_STAT
## AIC => 87.26753
##
## Table: Removing Existing Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## OUT_STAT 1 88.405 99.925 -121.264 0.44206 0.41815
## BOOK 1 89.065 100.585 -120.677 0.43706 0.41294
## PH_D 1 96.921 108.442 -113.661 0.37401 0.34718
## GRADUAT 1 109.030 120.550 -102.753 0.26272 0.23112
## -----------------------------------------------------------------------
##
##
## No more variables to be removed.
##
## Variables Removed:
##
## => log_ADD_FEE
## => SF_RATIO
## => R_B_COST
backward_aic
##
##
## Stepwise Summary
## --------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## --------------------------------------------------------------------------
## 0 Full Model 90.995 111.732 -117.098 0.48139 0.42639
## 1 log_ADD_FEE 89.001 107.434 -119.561 0.48135 0.43490
## 2 SF_RATIO 87.604 103.733 -121.356 0.47711 0.43866
## 3 R_B_COST 87.268 101.092 -122.021 0.46522 0.43422
## --------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ---------------------------------------------------------------
## R 0.682 RMSE 0.402
## R-Squared 0.465 MSE 0.162
## Adj. R-Squared 0.434 Coef. Var 12.482
## Pred R-Squared 0.385 AIC 87.268
## MAE 0.300 SBC 101.092
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 10.423 4 2.606 15.006 0.0000
## Residual 11.981 69 0.174
## Total 22.403 73
## -------------------------------------------------------------------
##
## Parameter Estimates
## ---------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ---------------------------------------------------------------------------------------
## (Intercept) -0.435 0.670 -0.649 0.519 -1.771 0.902
## BOOK 0.001 0.001 0.188 1.906 0.061 0.000 0.002
## PH_D 0.028 0.008 0.329 3.431 0.001 0.012 0.044
## GRADUAT 0.019 0.004 0.538 5.112 0.000 0.012 0.027
## OUT_STAT 0.000 0.000 -0.187 -1.729 0.088 0.000 0.000
## ---------------------------------------------------------------------------------------
forward_aic<-ols_step_forward_aic(model_pub_full, progress=TRUE, details=TRUE)
## Forward Selection Method
## ------------------------
##
## Candidate Terms:
##
## 1. BOOK
## 2. log_ADD_FEE
## 3. R_B_COST
## 4. PH_D
## 5. SF_RATIO
## 6. GRADUAT
## 7. OUT_STAT
##
##
## Step => 0
## Model => log_NEW10 ~ 1
## AIC => 125.5841
##
## Initiating stepwise selection...
##
## Table: Adding New Variables
## --------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## --------------------------------------------------------------------------
## GRADUAT 1 107.184 114.097 -103.864 0.24094 0.23040
## PH_D 1 107.752 114.664 -103.325 0.23509 0.22447
## R_B_COST 1 119.631 126.543 -92.031 0.10191 0.08943
## BOOK 1 121.310 128.222 -90.431 0.08130 0.06854
## log_ADD_FEE 1 122.064 128.977 -89.711 0.07188 0.05899
## SF_RATIO 1 126.507 133.419 -85.472 0.01445 0.00076
## OUT_STAT 1 127.284 134.197 -84.730 0.00404 -0.00979
## --------------------------------------------------------------------------
##
## Step => 1
## Added => GRADUAT
## Model => log_NEW10 ~ GRADUAT
## AIC => 107.1844
##
## Table: Adding New Variables
## -------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## PH_D 1 92.767 101.983 -117.502 0.39197 0.37484
## BOOK 1 96.918 106.134 -113.678 0.35689 0.33877
## OUT_STAT 1 104.070 113.286 -107.072 0.29163 0.27168
## R_B_COST 1 105.322 114.539 -105.913 0.27954 0.25924
## SF_RATIO 1 107.925 117.141 -103.501 0.25375 0.23273
## log_ADD_FEE 1 108.311 117.527 -103.142 0.24984 0.22871
## -------------------------------------------------------------------------
##
## Step => 2
## Added => PH_D
## Model => log_NEW10 ~ GRADUAT + PH_D
## AIC => 92.76671
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## BOOK 1 88.405 99.925 -121.264 0.44206 0.41815
## OUT_STAT 1 89.065 100.585 -120.677 0.43706 0.41294
## R_B_COST 1 94.202 105.723 -116.094 0.39659 0.37073
## log_ADD_FEE 1 94.654 106.174 -115.690 0.39290 0.36688
## SF_RATIO 1 94.759 106.280 -115.596 0.39203 0.36597
## ------------------------------------------------------------------------
##
## Step => 3
## Added => BOOK
## Model => log_NEW10 ~ GRADUAT + PH_D + BOOK
## AIC => 88.40477
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## OUT_STAT 1 87.268 101.092 -121.878 0.46522 0.43422
## R_B_COST 1 90.007 103.831 -119.523 0.44505 0.41288
## log_ADD_FEE 1 90.016 103.840 -119.515 0.44498 0.41281
## SF_RATIO 1 90.387 104.212 -119.195 0.44219 0.40986
## ------------------------------------------------------------------------
##
## Step => 4
## Added => OUT_STAT
## Model => log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT
## AIC => 87.26753
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## R_B_COST 1 87.604 103.733 -121.080 0.47711 0.43866
## SF_RATIO 1 88.775 104.904 -120.111 0.46877 0.42970
## log_ADD_FEE 1 89.109 105.237 -119.834 0.46636 0.42713
## ------------------------------------------------------------------------
##
##
## No more variables to be added.
##
## Variables Selected:
##
## => GRADUAT
## => PH_D
## => BOOK
## => OUT_STAT
forward_aic
##
##
## Stepwise Summary
## --------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## --------------------------------------------------------------------------
## 0 Base Model 125.584 130.192 -85.606 0.00000 0.00000
## 1 GRADUAT 107.184 114.097 -103.864 0.24094 0.23040
## 2 PH_D 92.767 101.983 -117.502 0.39197 0.37484
## 3 BOOK 88.405 99.925 -121.264 0.44206 0.41815
## 4 OUT_STAT 87.268 101.092 -121.878 0.46522 0.43422
## --------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ---------------------------------------------------------------
## R 0.682 RMSE 0.402
## R-Squared 0.465 MSE 0.162
## Adj. R-Squared 0.434 Coef. Var 12.482
## Pred R-Squared 0.385 AIC 87.268
## MAE 0.300 SBC 101.092
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 10.423 4 2.606 15.006 0.0000
## Residual 11.981 69 0.174
## Total 22.403 73
## -------------------------------------------------------------------
##
## Parameter Estimates
## ---------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ---------------------------------------------------------------------------------------
## (Intercept) -0.435 0.670 -0.649 0.519 -1.771 0.902
## GRADUAT 0.019 0.004 0.538 5.112 0.000 0.012 0.027
## PH_D 0.028 0.008 0.329 3.431 0.001 0.012 0.044
## BOOK 0.001 0.001 0.188 1.906 0.061 0.000 0.002
## OUT_STAT 0.000 0.000 -0.187 -1.729 0.088 0.000 0.000
## ---------------------------------------------------------------------------------------
Комментарий: строим Backward и Forward AIC. Backward AIC исключил всё, кроме SF_RATIO и GRADUAT. Forward AIC включил только GRADUAT и SF_RATIO. Оба пошаговых алгоритма дали один и тот же результат.
Интерпретация параметров: число лучших 10% студентов NEW10 увеличивается с увеличением выпускников GRADUAT, числом докторов (PH_D), ценой за бронирование (BOOK) и уменьшается со стоимостью оплаты OUT_STAT.
Здесь можно сравнить по R^2, например.
model_pub_full<-lm(log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT, data=pub_data)
summary(model_pub_full)
##
## Call:
## lm(formula = log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.21130 -0.24350 -0.01198 0.24097 0.92728
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.967e-02 5.542e-01 -0.108 0.91449
## GRADUAT 1.825e-02 3.483e-03 5.241 9.95e-07 ***
## PH_D 2.271e-02 7.077e-03 3.209 0.00183 **
## BOOK 1.152e-03 4.545e-04 2.535 0.01290 *
## OUT_STAT -3.446e-05 2.284e-05 -1.509 0.13474
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4289 on 93 degrees of freedom
## (25 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4553, Adjusted R-squared: 0.4319
## F-statistic: 19.43 on 4 and 93 DF, p-value: 1.196e-11
AIC(model_pub_full)
## [1] 119.0812
summary(model_pub_nono)
##
## Call:
## lm(formula = log_NEW10 ~ R_B_COST + PH_D + GRADUAT + OUT_STAT +
## BOOK, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.21996 -0.24955 -0.00251 0.24842 0.81900
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.390e-02 5.501e-01 0.134 0.8934
## R_B_COST 1.272e-04 6.475e-05 1.964 0.0525 .
## PH_D 1.748e-02 7.462e-03 2.342 0.0213 *
## GRADUAT 1.808e-02 3.432e-03 5.268 9.03e-07 ***
## OUT_STAT -5.027e-05 2.389e-05 -2.104 0.0381 *
## BOOK 1.044e-03 4.510e-04 2.315 0.0228 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4225 on 92 degrees of freedom
## (25 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4772, Adjusted R-squared: 0.4488
## F-statistic: 16.8 on 5 and 92 DF, p-value: 9.211e-12
AIC(model_pub_nono)
## [1] 117.0556
Комментарий: эта модель в сравнении с model_pub_nono имеет более высокий p-value и более низкий R^2.
plot(model_pub_full, 1:5)
plot(residuals(model_pub_full), rstudent(model_pub_full))
Комментарий: особых выбросов нет.
influence.measures(model_pub_full)
## Influence measures of
## lm(formula = log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT, data = pub_data) :
##
## dfb.1_ dfb.GRAD dfb.PH_D dfb.BOOK dfb.OUT_ dffit cov.r cook.d
## 1 -0.042131 -0.133457 3.02e-02 0.016099 1.27e-01 -0.17853 1.047 6.39e-03
## 3 0.103254 0.106100 -1.30e-01 -0.039284 4.56e-02 -0.20747 1.147 8.67e-03
## 5 0.025264 0.031849 -1.62e-02 -0.037843 -1.18e-02 -0.07080 1.072 1.01e-03
## 6 0.042930 -0.026515 -3.75e-02 0.020769 2.19e-03 0.07699 1.074 1.20e-03
## 7 -0.001421 -0.001482 2.05e-03 -0.000293 3.65e-05 0.00295 1.083 1.76e-06
## 8 0.256242 -0.005185 -1.43e-01 -0.095065 -1.30e-01 0.30174 1.034 1.81e-02
## 9 -0.223857 0.071406 9.59e-02 0.076188 2.17e-01 0.40619 0.932 3.22e-02
## 11 -0.350149 -0.136018 1.30e-01 0.312988 3.54e-01 0.52495 1.037 5.43e-02
## 12 -0.121942 0.075940 3.86e-02 0.138385 -2.94e-02 0.21075 1.137 8.94e-03
## 13 -0.215844 -0.020160 1.53e-01 0.134773 -1.49e-03 0.27751 1.055 1.54e-02
## 15 -0.321287 -0.077636 2.49e-01 0.025996 2.44e-01 0.44045 0.919 3.77e-02
## 16 -0.250015 -0.049278 2.92e-01 -0.032726 -4.84e-02 0.38214 0.833 2.80e-02
## 17 -0.081306 0.044052 6.07e-02 0.007316 -2.62e-02 -0.11473 1.077 2.65e-03
## 18 0.003194 -0.000841 1.97e-02 -0.041553 -4.86e-03 0.05185 1.085 5.43e-04
## 19 0.009679 0.007517 -7.39e-03 0.003912 -2.18e-02 -0.02766 1.122 1.55e-04
## 20 0.141867 0.108893 -1.32e-01 0.020169 -1.07e-01 0.23226 0.956 1.06e-02
## 21 -0.012381 0.071752 1.14e-02 0.006585 -6.69e-02 0.12529 1.036 3.15e-03
## 22 0.036851 -0.032196 3.27e-02 -0.086552 -3.76e-02 0.12703 1.053 3.24e-03
## 23 -0.112694 0.134024 1.63e-02 0.188038 -1.08e-01 0.30489 1.045 1.85e-02
## 24 0.134601 0.381840 -1.22e-01 -0.241573 -1.15e-01 -0.53441 0.843 5.46e-02
## 25 0.014454 0.059799 3.86e-02 -0.077258 -1.06e-01 0.13967 1.084 3.93e-03
## 26 0.006213 0.014783 -3.94e-03 0.011396 -3.38e-02 0.04477 1.100 4.05e-04
## 27 -0.003111 0.012714 4.84e-02 -0.119863 2.89e-03 -0.13920 1.130 3.91e-03
## 28 -0.224463 -0.046283 1.85e-01 0.043128 1.35e-02 -0.27959 0.924 1.53e-02
## 29 -0.199615 -0.064122 1.55e-01 0.074507 2.84e-02 -0.22725 1.057 1.03e-02
## 31 0.024836 0.118218 9.65e-05 -0.070988 -8.35e-02 0.15237 1.108 4.68e-03
## 32 0.010801 -0.068147 -4.39e-04 0.018210 2.63e-02 0.08799 1.083 1.56e-03
## 33 0.546814 0.131824 -4.82e-01 -0.062290 -6.47e-02 0.58654 0.888 6.62e-02
## 34 -0.017402 -0.014653 2.72e-02 -0.010574 -9.31e-03 -0.03488 1.106 2.46e-04
## 35 0.000434 -0.007105 1.32e-03 0.002973 -6.78e-03 -0.02361 1.072 1.13e-04
## 36 -0.030540 -0.044774 4.61e-02 -0.026909 1.86e-02 -0.07377 1.075 1.10e-03
## 38 0.038315 0.014707 -3.98e-02 0.010623 -9.85e-03 0.05040 1.086 5.13e-04
## 39 0.031172 -0.009854 -1.04e-02 -0.062093 2.45e-02 -0.11232 1.044 2.53e-03
## 41 -0.001574 0.010763 -4.89e-03 0.005513 -4.69e-04 -0.01433 1.110 4.15e-05
## 42 0.033069 -0.087660 1.88e-02 -0.011504 -2.87e-02 0.15290 1.040 4.69e-03
## 44 0.043104 0.028398 -9.65e-02 0.086072 1.48e-02 -0.12728 1.088 3.26e-03
## 45 0.027188 -0.009637 -5.15e-02 0.053262 -1.60e-03 -0.11268 1.036 2.55e-03
## 46 0.032962 -0.006529 -6.83e-03 -0.044942 8.94e-03 0.07515 1.064 1.14e-03
## 47 -0.130253 0.121583 6.48e-02 0.136264 -2.07e-01 -0.38168 0.845 2.80e-02
## 48 -0.000454 -0.146705 -1.18e-01 0.270710 8.78e-02 -0.40760 0.809 3.17e-02
## 49 0.091508 -0.006039 -1.23e-01 0.115218 -6.34e-02 -0.22342 1.046 9.98e-03
## 51 0.020780 0.007590 -8.55e-03 -0.017144 -1.30e-02 0.03116 1.082 1.96e-04
## 52 0.005839 -0.001717 -4.46e-03 -0.003079 4.69e-03 0.01013 1.100 2.07e-05
## 53 -0.002512 -0.001853 1.40e-03 0.000344 3.86e-03 -0.00499 1.094 5.04e-06
## 54 0.017187 0.003092 -5.29e-03 -0.012185 -1.72e-02 0.02580 1.111 1.35e-04
## 55 -0.007263 -0.006925 6.30e-03 0.000973 1.36e-02 0.03071 1.066 1.91e-04
## 56 -0.102366 -0.142662 4.14e-02 0.142251 1.49e-01 0.25798 0.993 1.32e-02
## 58 0.071336 0.007929 -5.81e-02 -0.012835 -6.21e-03 0.07724 1.110 1.20e-03
## 59 -0.020325 0.005678 1.20e-02 0.003531 7.44e-03 -0.02635 1.101 1.40e-04
## 61 0.026528 0.015272 3.82e-02 -0.084232 -8.32e-02 0.12538 1.081 3.17e-03
## 63 0.087394 -0.059968 1.25e-02 -0.090035 -1.38e-01 -0.24911 1.013 1.23e-02
## 64 -0.170093 -0.191554 2.91e-01 -0.103360 -5.94e-02 0.34442 1.146 2.38e-02
## 66 0.247524 -0.418324 -1.60e-01 -0.152779 3.84e-01 -0.69659 0.693 8.92e-02
## 67 0.034292 0.137055 -7.46e-02 0.059278 -8.08e-02 0.16454 1.130 5.46e-03
## 68 0.026783 -0.003626 -5.07e-02 0.024273 4.16e-02 -0.08292 1.068 1.39e-03
## 69 -0.019209 -0.001041 2.31e-02 -0.007225 -4.35e-04 0.03438 1.075 2.39e-04
## 70 0.012013 0.239901 -9.68e-03 -0.078056 -1.17e-01 0.28916 1.023 1.66e-02
## 71 -0.092085 -0.005378 1.14e-01 -0.044926 -4.20e-02 -0.14123 1.070 4.01e-03
## 72 0.093203 0.066536 -2.64e-01 0.242285 1.34e-01 -0.34264 1.152 2.35e-02
## 73 -0.158831 -0.232355 1.86e-01 -0.033745 1.25e-01 -0.32859 0.894 2.10e-02
## 74 -0.110318 0.092402 9.43e-03 0.127219 -2.19e-02 -0.23596 0.963 1.10e-02
## 75 -0.062227 -0.055431 5.68e-02 0.012618 2.93e-02 -0.09321 1.070 1.75e-03
## 77 0.002545 0.003570 -4.48e-03 0.002996 -2.57e-03 -0.00797 1.079 1.28e-05
## 79 0.000397 -0.000864 1.38e-02 -0.043803 1.39e-02 -0.05684 1.119 6.53e-04
## 80 -0.003717 -0.005744 -5.30e-04 0.022151 -1.00e-02 0.03494 1.106 2.47e-04
## 83 0.030763 -0.027421 -4.54e-02 0.016224 7.42e-02 0.10076 1.086 2.05e-03
## 84 0.051012 0.143654 -5.12e-02 -0.013781 -1.43e-01 -0.19210 1.059 7.40e-03
## 85 0.023201 0.012497 -5.37e-02 0.070858 -2.02e-02 -0.09385 1.147 1.78e-03
## 86 -0.045718 0.003991 1.31e-02 0.126471 -1.41e-01 -0.29579 0.921 1.71e-02
## 87 -0.009540 -0.062359 6.97e-02 -0.098075 2.55e-03 -0.12517 1.147 3.16e-03
## 88 0.014200 0.014410 -2.16e-04 -0.027171 -1.11e-02 0.04109 1.078 3.41e-04
## 90 0.008670 -0.012445 1.30e-02 -0.071909 3.83e-02 -0.10831 1.080 2.36e-03
## 91 0.029025 -0.048617 7.56e-02 -0.120907 -9.13e-02 0.17578 1.144 6.23e-03
## 92 0.430023 -0.528612 -4.13e-01 0.206695 3.71e-01 0.83959 1.084 1.38e-01
## 93 0.045189 -0.017476 1.07e-02 -0.050952 -6.48e-02 0.10855 1.083 2.37e-03
## 96 0.196067 -0.069055 -8.44e-02 -0.082115 -9.60e-02 0.26129 1.081 1.37e-02
## 97 -0.009366 0.027867 1.62e-02 -0.002606 -4.52e-02 0.05736 1.107 6.65e-04
## 98 0.027003 0.036114 -7.06e-03 -0.057481 -1.77e-02 -0.07484 1.166 1.13e-03
## 100 0.075991 0.026069 1.59e-02 0.035124 -3.98e-01 -0.48472 1.107 4.67e-02
## 101 -0.028942 0.123097 -2.87e-02 0.017744 3.96e-02 0.18388 1.111 6.80e-03
## 102 -0.004197 0.007296 5.48e-03 -0.003717 -9.50e-03 -0.01427 1.096 4.12e-05
## 103 0.007016 -0.020480 -1.38e-03 -0.009095 2.28e-02 0.03411 1.105 2.35e-04
## 104 0.123240 0.179060 -1.41e-01 0.026554 -1.63e-01 -0.24651 1.098 1.22e-02
## 105 0.023884 -0.063654 8.23e-02 -0.146901 -7.63e-02 -0.20001 1.083 8.03e-03
## 106 0.055488 -0.001506 -4.34e-02 -0.022109 5.94e-04 -0.07177 1.090 1.04e-03
## 107 0.004810 0.000641 1.61e-02 -0.041155 -1.16e-02 -0.04542 1.121 4.17e-04
## 109 0.015535 -0.011063 -2.23e-02 0.019857 3.92e-03 -0.03938 1.093 3.13e-04
## 110 -0.045418 0.099942 3.16e-02 -0.007424 -8.10e-02 -0.14449 1.071 4.20e-03
## 111 0.003248 0.005327 -1.45e-03 -0.004667 -4.85e-03 -0.00846 1.105 1.45e-05
## 113 0.037474 -0.129469 3.64e-03 -0.054592 7.57e-02 -0.17693 1.065 6.28e-03
## 114 0.223708 0.055697 -2.56e-01 -0.035107 2.07e-01 0.42656 0.901 3.53e-02
## 115 -0.039680 0.195298 2.60e-03 -0.087338 5.77e-02 0.32221 1.088 2.07e-02
## 116 -0.012937 -0.026990 1.32e-02 0.008373 9.61e-03 -0.03206 1.140 2.08e-04
## 117 -0.137025 0.039881 4.03e-02 -0.106701 4.40e-01 0.61275 1.075 7.40e-02
## 118 0.001980 -0.000154 -1.45e-03 -0.002237 1.56e-03 -0.00497 1.084 4.99e-06
## 119 -0.195111 0.041427 1.92e-01 -0.064464 -1.77e-02 -0.25500 1.043 1.30e-02
## 122 -0.049622 -0.115513 9.73e-02 -0.034789 2.04e-02 0.14702 1.104 4.35e-03
## 123 0.002463 0.058149 1.39e-02 -0.020111 -7.73e-02 0.09267 1.112 1.73e-03
## hat inf
## 1 0.0367
## 3 0.0988
## 5 0.0252
## 6 0.0278
## 7 0.0254
## 8 0.0578
## 9 0.0482
## 11 0.1039
## 12 0.0931
## 13 0.0613
## 15 0.0513
## 16 0.0283 *
## 17 0.0378
## 18 0.0312
## 19 0.0595
## 20 0.0232
## 21 0.0210
## 22 0.0278
## 23 0.0626
## 24 0.0515
## 25 0.0471
## 26 0.0428
## 27 0.0772
## 28 0.0256
## 29 0.0513
## 31 0.0647
## 32 0.0365
## 33 0.0693
## 34 0.0470
## 35 0.0171
## 36 0.0278
## 38 0.0325
## 39 0.0210
## 41 0.0491
## 42 0.0283
## 44 0.0472
## 45 0.0181
## 46 0.0214
## 47 0.0296
## 48 0.0288 *
## 49 0.0460
## 51 0.0267
## 52 0.0403
## 53 0.0354
## 54 0.0508
## 55 0.0137
## 56 0.0355
## 58 0.0543
## 59 0.0418
## 61 0.0420
## 63 0.0393
## 64 0.1210
## 66 0.0505 *
## 67 0.0814
## 68 0.0257
## 69 0.0209
## 70 0.0510
## 71 0.0391
## 72 0.1239
## 73 0.0286
## 74 0.0248
## 75 0.0290
## 77 0.0218
## 79 0.0593
## 80 0.0472
## 83 0.0403
## 84 0.0447
## 85 0.0845
## 86 0.0277
## 87 0.0877
## 88 0.0242
## 90 0.0379
## 91 0.0922
## 92 0.1798 *
## 93 0.0399
## 96 0.0703
## 97 0.0495
## 98 0.0974 *
## 100 0.1264
## 101 0.0723
## 102 0.0374
## 103 0.0464
## 104 0.0760
## 105 0.0584
## 106 0.0386
## 107 0.0597
## 109 0.0362
## 110 0.0403
## 111 0.0452
## 113 0.0442
## 114 0.0452
## 115 0.0860
## 116 0.0745
## 117 0.1361
## 118 0.0261
## 119 0.0515
## 122 0.0610
## 123 0.0579
# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pub_full)
cook_df <- data.frame(
obs = row.names(model_pub_full$model),
cook = cooks_vals
) %>%
arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))
plot(
x, cook_df$cook,
type = "h",
main = "Cook's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Cook's distance",
xaxt = "n" # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)
Комментарий: высокие значения у 29 и 38.
# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pub_full)[,-1],
center = colMeans(model.matrix(model_pub_full)[,-1]),
cov = cov(model.matrix(model_pub_full)[,-1]))
mah_df <- data.frame(
obs = row.names(model_pub_full$model),
mah = mahalanobis_distance
) %>%
arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05
# Scatter plot для расстояния Махаланобиса
plot(
x, mah_df$mah,
type = "h",
main = "Mahalanobis's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Mahalanobis's distance",
xaxt = "n"
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)
set.seed(456) # Другое зерно для генерации новых значений
max_GRADUAT<-max(pub_data$GRADUAT, na.rm=TRUE)
max_BOOK<-max(pub_data$BOOK, na.rm=TRUE)
max_OUT_STAT<-max(pub_data$OUT_STAT, na.rm=TRUE)
max_PH_D<-max(pub_data$PH_D, na.rm=TRUE)
new_data <- data.frame(
GRADUAT = max_GRADUAT - 1:8,
BOOK = max_BOOK - 1:8,
OUT_STAT = max_OUT_STAT - seq(100, 800, by=100),
PH_D = max_PH_D - 1:8
)
# Прогнозируем значение NEW10
predicted_values <- predict(model_pub_full, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")
## Прогнозируемые значения NEW10:
print(predicted_values)
## 1 2 3 4 5 6 7 8
## 4.330459 4.291791 4.253122 4.214454 4.175785 4.137117 4.098448 4.059780
# 6. Доверительные и предсказательные интервалы
# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pub_full, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")
## [1] "Доверительные интервалы (95%):\n"
print(confidence_intervals)
## fit lwr upr
## 1 4.330459 3.917940 4.742979
## 2 4.291791 3.883437 4.700145
## 3 4.253122 3.848577 4.657668
## 4 4.214454 3.813349 4.615558
## 5 4.175785 3.777744 4.573826
## 6 4.137117 3.741753 4.532480
## 7 4.098448 3.705369 4.491527
## 8 4.059780 3.668584 4.450975
# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pub_full, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")
## [1] "Предсказательные интервалы (95%):\n"
print(prediction_intervals)
## fit lwr upr
## 1 4.330459 3.384023 5.276896
## 2 4.291791 3.347162 5.236419
## 3 4.253122 3.310134 5.196111
## 4 4.214454 3.272936 5.155971
## 5 4.175785 3.235569 5.116001
## 6 4.137117 3.198031 5.076202
## 7 4.098448 3.160322 5.036574
## 8 4.059780 3.122441 4.997118
# Создаем data.frame для прогнозов и интервалов
predictions_df <- data.frame(
GRADUAT = new_data$GRADUAT,
OUT_STAT = new_data$OUT_STAT,
BOOK = new_data$BOOK,
PH_D = new_data$PH_D,
Predicted = predicted_values,
Lower_Conf = confidence_intervals[, "lwr"],
Upper_Conf = confidence_intervals[, "upr"],
Lower_Pred = prediction_intervals[, "lwr"],
Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df
# Визуализация (только для демонстрации, нужен более сложный график для нескольких предикторов)
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) + # Условная ось x, т.к. несколько предикторов
geom_point() +
geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
x = "Прогноз (номер)",
y = "NEW10") +
scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
theme_bw()
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.9, label = "Синие: Доверительные интервалы", color = "blue", hjust = 0) +
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.8, label = "Красные: Предсказательные интервалы", color = "red", hjust = 0)